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ABSTRACT 


Two  algorithms  for  generating  a  non-homogeneous 
Poisscn  process  with  log-quadratic  intensity  function 
A  (t)  =  exp  (cxq  ♦  ait  ♦  a2"t  )  are  implemented  into 
computer  programs  and  compared  for  relative  speed, 
core  storage  requirements  and  fidelity.  3y  simulating 
several  cases  of  non-homogeneous  Poisscn  processes 
with  log-quadratic  intensity  functions  it  is  shewn 
that  -che  Poisson-decomposition  and  gap  statistic 
algorithm,  developed  by  Professor  P.A.W.  Lewis,  Naval 
Postgraduate  School,  Monterey,  California,  and  G.S. 
Shedler,  IBM  Research  Laboratory,  San  Jose, 
California,  substantially  reduces  computation  time 
from  that  required  by  an  algorithm  that  uses  a 
time-scale  transformation  of  a  homogeneous  Poisson 
process.  The  faster  algorithm  employs  a  rejection 
technigue  in  conjunction  with  a  method  for  simulating 
the  ncn-hc mogeneous  Poisson  process  with  intensity 
function  A(t)  =  exp  (  y  +  y  t)  by  generation  of  gap 
statistics.  Although  additional  core  storage  is 
required  by  the  Lewis  and  Shedler  algorithm,  the 
resulting  gain  in  computing  efficiency  is  so 
significant  that  it  outweighs  the  memory 
consideration.  The  experience  gained  from 
implementing  the  algorithm  has  led  to  several 
possibilities  which  are  suggested  for  improving  the 
efficiency  of  the  Poisson-decomposition  and  gap 
statistic  algorithm. 
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I.   INTRODUCTION 


Many  familiar  physical  and  operational  processes  are 
well  described,  in  whole  or  part,  by  examining  their  "event 
streams"  ever  time.  Some  such  well-known  processes  are:  the 
flow  of  traffic  through  an  intersection;  the  arrival  of 
persons  at  seme  service  facility  such  as  a  bank  teller's 
window,  a  service  station  fuel  pump  or  a  grocery  store  check 
out  counter;  and,  the  arrival  of  telephone  calls  or  radio 
transmissions  at  some  switchboard  or  other  type  of 
communications  terminal.  Analogous  processes  abound  both  in 
nature  and  in  the  course  of  our  everyday  lives.  Ey  the 
proper  definition  of  an  "event"  in  these  situations,  the 
process  being  observed  will  be  characterized  by  the 
probabilistic  nature  of  the  flow,  over  time,  of  the  events 
of  which  it  is  composed.  In  the  above  three  examples  the 
events  could  be  defined  respectively  as  the  arrival  of  a 
vehicle  at  the  intersection,  the  arrival  of  a  customer  at 
the  service  facility,  and  the  arrival  of  the  telephone  call 
or  radio  transmission  at  the  terminal.  The  process  may  then 
be  analyzed  ty  examining  the  interaction  of  various  event 
streams  with  different  intersection  configurations,  service 
policies  and  terminal  capacities. 

A  common  method  used  to  perform  such  analyses  is  to 
consider  the  event  streams  to  be  homogeneous.  This  could 
mean  that  the  expected  number  of  events  to  occur  in  any  two 
or  more  time  intervals  of  equal  length  is  the  same  (simple 
homogeneity) ,  or  that  the  distribution  of  the  number  of 
events  occurring  in  any  two  or  more  equal  time  intervals  is 
the  same  (complete  homogeneity) .  The  homogeneous  Poisson 
process  is  often  used  as  a  tool   in   the   analysis   of   such 
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activities.  For  very  simple  systems  involving  event  streams 
the  use  of  the  homogeneous  Pcisson  process  as  a  model  leads 
to  tractable  analytical  results.  Most  systems  of  interest, 
however,  are  not  amenable  to  purely  analytical  methods,  and 
simulation  of  the  processes  by  digital  means  becomes 
necessary. 

The  assumption  of  homogeneity  in  event  streams  is  often 
a  very  restrictive  one.  The  "rush  hour"  phenoienon, 
well-known  and  abundantly  cursed  by  motorists,  provides 
cogent  evidence  that  the  modeling  of  event  streams  is  not 
always  well  served  by  the  homogeneity  assumption.  The 
intensity  of  event  streams  varies  over  time  for  many 
physical  cr  operational  processes.  In  these  cases,  purely 
analytical  methods  must  be  abandoned  almost  immediately  in 
fawcr  of  simulation  techniques.  (The  intensity  of  the  event 
stream  is  defined  to  be  the  derivative  with  respect  to  t  of 
the  expected  number  of  events  in  an  interval  of  length 
(0,t].) 

The  ncn-hcmogeneous  Poisscn  process  is  often  employed  in 
the  analysis  of  processes  that  exhibit  gross  departures  from 
the  homogeneous  event  stream  criterion.  If  these  processes 
are  to  be  simulated,  it  is  necessary  to  first  describe  the 
nature  of  the  inhoraogeneity  (i.e.  how  does  the  intensity  of 
the  event  stream  vary  over  time?)  and  then  to  artificially 
generate  event  streams  that  behave  in  accordance  with  the 
descripticn. 

Of  course  there  are  infinite  variations  in  the  types  of 
inhomogeneities  that  can  occur  or  that  can  be  construed. 
However  it  is  intuitively  appealing  to  consider  event 
streams  that  display  varying  intensities  of  the  following 
types: 

i)  increasing  continuously  over  time; 
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ii)  decreasing  continuously  over  time; 

iii)  increasing  and   then   decreasing   continuously   over 
time ; 

iv)  decreasing  then  increasing  continuously  over  time. 

A  ncn-hcmogeneous  Poisson  process  that  has  guadratic 
properties  can  be  manipulated  to  produce  the  above-mentioned 
effects.  The  effective  simulation  of  such  a  process  on  the 
computer  is  the  subject  of  this  thesis  and  has  been 
motivated  by  the  work  of  P.  A.  W.  Lewis,  Professor,  Naval 
Postgraduate  School,  and  G.  S.  Shedler,  IBM  Research 
Laboratory. 

There  are  of  course  other  types  of  inhomogeneities ,  such 
as  cyclic  variations  (time  of  day  effect) ,  but  we  do  not 
consider  them  here.  They  have  been  discussed  by  COX 
[Sef.  2]  and  LEWIS  [Ref.  9].  In  particular  LEWIS  [Ref.  9] 
describes  a  process  consisting  of  arrivals  at  an  intensive 
care  unit  in  a  hospital.  It  is  shown  empirically  that,  in 
addition  to  the  time  of  day  cycle,  long  term  fluctuations  in 
the  intensity  function  can  be  adequately  described  by  an 
intensity  function  whose  logarithm  is  quadratic. 

LEWIS  and  SHEDLER  [Ref.  11]  proposed  a  new  method  for 
generating  a  non-homogeneous  Poisson  process  with  an  event 
stream  intensity  (rate)  function  that  is  of  degree-two 
exponential  polynomial  form.  (The  use  of  exponential 
polynomials  is  natural  in  this  context  since  an  intensity 
function  is  a  positive  function.)  The  new  method  appears  to 
have  the  virtue  of  increased  efficiency  over  the  more 
conventional  time-scale  transformation  technique  when 
implemented  on  a  high  speed  digital  computer. 

Efficiency  in  this   context   is   measured   in   terms   of 
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computer  memory  requirements  and  computational  speed.  The 
problem  of  efficiency  comparison  is  recognized  by  LEWIS  and 
SHZDLER  in  the  final  pages  of  their  paper  [ Ref .  11 f  p.  15]: 
"There  remains  the  question  of  efficiency  (cf  the  proposed 
algorithm)  .  .  .  for  generation  of  a  non-homogeneous  Pcisson 
process  with  degree-two  exponential  polynomial  rate 
function,  relative  to  generation  via  time  scale 
transformation  of  a  homogeneous  Poisson  process." 

After  seme  brief  discussion  of  the  requirements  of 
inplementing  the  time-scale  transformation  algorithm  the 
report  concludes  .  .  .  "We  therefore  would  expect  the  exact 
method  of  (the  proposed  algorithm)  to  be  much  faster, 
although   at  the  expense  of  some  complexity  of  programming." 

The  objective  of  this  author  has  been:  to  implement 
both  the  algorithm  of  LEWIS  and  SHEDLER  [Ref.  11]  and  the 
conventional  time-scale  transformation  algorithm  on  the  IBM 
360/67  Computer  System  in  FORTRAN  IV  language;  to  define 
reasonable  measures  of  effectiveness  for  comparing  th€  two 
algorithms;  and  to  determine  which  algorithm  is  the  more 
efficient  in  terras  of  the  measures  of  effectiveness  defined. 

Section  II  discusses  the  definition  of  a  non-homogeneous 
Pcisson  process.  It  also  states  some  special  properties  of 
Pcisson  processes  that  are  used  in  the  development  cf  the 
algorithms  investigated  in  this  thesis,  and  concludes  sith  a 
general  discussion  of  two  basic  methods  of  generating  a 
non-homogenecus  Poisson  process.  Section  III  gives  a  step 
by  step  description  of  the  two  algorithms  that  were 
implemented  into  computer  programs.  The  method  of 
implementation  is  the  topic  of  Section  IV.  Methods  of 
comparing  the  algorithms  are  presented  in  Section  V. 
Section  VI  presents  conclusions  drawn  from  the  comparison 
and  makes  a  recommendation  for  improving  the  LEWIS  and 
SHEDLER   [Ref.   11]   algorithm.    Other   recommendations  for 
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further  study  are  also  listed  in  Section  VI.  Appendixes 
provide  additional  details  that  would  have  been  awkward  to 
include  in  the  main  body  of  the  thesis.  Computer  program 
listings  are  provided  after  Appendix  E. 
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II.        DEFINITION    AND    PROPERTIES    OF    NON- HOMOGENEOUS    PCISSON 

PROCESSES 


A.       GENERAL 


This  section  will  present  basic  definitions  and 
explanations  concerning  the  concept  of  non-homogeneous 
Pcisson  processes.  References  cited  may  be  consulted  if  a 
mere      in-depth        understanding        is        desired.  Only        the 

fundamental  concepts  necessary  for  understanding  the 
specific  non-homogeneous  Poisson  process  under  consideration 
are    presented. 


B.       DEFINITION    OF    A    NON-HOMOGENEOUS    POISSON    PROCESS 


LEWIS    and   SHEDLER    [Ref.     11]    define      the      non-homogeneous 
Poisson    process   on    the   real    line    as   follows: 

1.  The      number      of      events      in      any         finite         set        of 
non-overlapping      intervals    are    independent    random    variables. 

2.  Let        A (t)         be  a  monotone  non- deer  easing 

right-continuous      function      which      is      bounded   on    any    finite 

interval.      Then   the    number    of    events    in    any      interval,       e.g. 

(0,tQ),         has        a         Poisson         distribution         with       parameter 

A(tQ)     -      A  (0)  . 

The    function      A (t)     -   A (0)     is    called,    among    other   things,    the 
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mean  value  function  since  the  expected  number  of  events  in 
an  interval  (0,t]  is  just  A  (t)  -  A  (0)  .  Property  1  is  the 
independent  increment  property  of  the  Poisson  process;  it  is 
basic  to  the  idea  of  a  Poisson  process.  Figure  1  provides  a 
graphical  representation  of  the  definition  above. 

The       above       definition       insures       that 

(A(t)  -   A(0)}/{  A(t  )  -   A(0)}    meets     the     following 

criteria   for   an   arbitrary   function   F  (t) ,   to  be  a  valid 

distribution  function  on  the  interval  (0,tQ),   c.f.    LARSON 

[Ref.  7,  ch.  3]; 

i)   0  <_  F(t)  <_   1 

ii)    lim  F(t)  =  0;   lim   F(t)  -   1,       0  <  t  <  t 
t  +  0  t  -  tQ 

iii)   F(a)  <_  F(b)   for  all   a  <_   b   in   (0,t  ) 

iv)   lim  F(b  +  h)  =  F(b)   for  all   b   in   (0,tn), 
h  +  0  u 

where   h  >  0  . 


Letting  F  (t)  =  {  A  (t)  -  A(0)}  /  C  A  (t  )  -  A  (0)  }  it  follows 
that  if  A (t)  is  absolutely  continuous  in  (0,t  ]  then 
dF(t)/dt  =  A(t)/{  A  (t  )  -  A(0)}  is  a  valid  density  function 
on  the  interval  (0,t  ].  The  function  \  (t)  is  called  the 
intensity  function  (or  rate  function)  of  the  process  and 


X (t)  =  ^-  E{number  of  events  in   (0,tfl]} 

=  ^  CA(t)  -  A(0)} 

d 
=  dt  A(t)  . 
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♦    A(t) 


Situation:      Events  occur  randomly   in   time   (i.e.  along   t-axis) 

1)  If     X  =  number  of  events   in  fixed   interval  (t, ,t_] 

Y   =   number  of  events    in   fixed   interval  (t?,t-J 

Z   =  number  of  events    in   fixed   interval  Uvt-] 

S  =  number  of  events   in  fixed   interval  (O.t-J; 

Then  X,  Y  and  Z  must  be  independent  random  variables  (note: 

S  and  Z  are  not  independent  because  of  overlap  in  the  intervals). 

2)  Given  the  monotonic  increasing  right  continuous  function 

A(t),  the  number  of  events  in  any  fixed  interval   (t-,t  1 

i      J 

must  have   the  Poisson  distribution  with  parameter     A(t.)    -  Aft.) 

J  l 

Thus 

X  ~  Poisson{A(t2)   -  A(t1)} 

Y  ~  Poisson(A(t3)   -  A(t2)} 

Z  ~  Poisson{A(t4)  -  A(t3)} 

S  -  Poisson{A(tQ)  -  a(0)>   . 


Fiaure    1    - 


DEFINITION    OF     A    NON-HOMOGENEOUS    POISSON    PROCESS 
(GRAPHICAL    REPRESENTATION) 
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THE  INTENSITY  FUNCTION 


The  most  intuitively  appealing  way  to  think  about  a 
non-homogeneous  Poisson  process  is  to  consider  the  variation 
in  the  intensity  of  the  event  stream  over  time.  The  concept 
of  an  intensity  function  whose  integral  meets  the  criteria 
of  the  definition  in  paragragh  B  above  is  essential  tc  the 
modeling  and  simulation  of  a  non-homogeneous  Poisson 
process.  When  one  starts  with  the  existence  of  A  (t) ,  the 
A(t)  is  often  called  the  integrated  intensity  (or  rate) 
function. 

The  intensity  functuion  reveals  the  instantaneous  rate 
of  arrivals  (in  the  event  stream)  as  a  function  of  time. 
(This  function  must  be  a  positive  function.)  For  example, 
if  telephone  calls  arrive  at  a  switchboard  at  the  rate  of  5 
per  hour  at  0900  (9:00  a.m.)  ,  increase  to  a  peak  rate  of  20 
per  hour  at  1300  (1:00  p.m.)  ,  then  decrease  to  a  rate  of  5 
per  hour  at  1700  (5:00  p.m.),  the  intensity  function  could 
lock  something  like  Figure  2a.  Then,  by  plotting  the 
integral  of  the  intensity  function  over  the  interval  of 
interest  (i.e.  from  0900  to  1700)  it  is  obvious  that  a 
monotone-increasing,  right-continuous  and  bounded  function 
is  obtained  (see  Figure  2b).  If  the  assumption  is  made  that 
the  arrival  stream  is  a  Poisson  process,  i.e.  has 
independent  increments,  then  the  number  of  calls  received  in 
any  chosen  interval  (e.g.  0900-1000,  1230-1315,  1107-1632) 
is  distributed  as  a  Poisson  random  variable  with  parameter 
egual  tc  the  difference  between  the  integral  evaluated  at 
the  right  end  point  of  the  interval  and  the  integral 
evaluated  at  the  left  end  point  of  the  interval.  These 
values  may  be  read  directly  from  Figure  2b.  Specifically, 
the  number  of  calls  received  in  an  eight-hour  working  day  is 
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a  Scissor  random  variable  with  parameter  120. 

Although  it  is  unlikely  that  telephone  calls  would 
arrive  at  a  switchboard  in  accordance  with  the  convenient 
parabolic  intensity  function  of  Figure  2a,  the  example 
serves  to  illustrate  two  important  points.  First,  it  would 
not  be  realistic  to  assume  that  the  arrival  rate  of 
telephone  calls  at  a  switchboard  would  be  constant 
throughout  a  working  day.  Some  function  that  describes  an 
initially  increasing  and  finally  decreasing  rate  of  arrivals 
seems  more  akin  to  reality.  The  importance  of  being  able  to 
model  a  non-homogeneous  Poisson  process  is  therefore 
established. 

Secondly,  it  is  obvious  that  the  definition  of  a 
ncn-homogeneous  Poisson  process  is  not  always  used  as  a 
starting  point  fcr  modeling  operational  processes.  Svent 
streams  are  usually  thought  of  in  terms  of  their  underlying 
intensity  functions.  The  idea  of  an  intensity  function 
applies  to  any  model  for  an  event  stream,  and  is  not 
specific  to  a  Pcisson  process.  The  further  step  of  modeling 
the  physical  process  as  a  non-homogeneous  Poisson  process  by 
assuming  that  the  process  has  independent  increments  is 
taken  either  on  the  basis  of  empirical  evidence  or  physical 
reasoning.  Testing  for  independent  increments  in  a  point 
process  is  discussed  in  COX  and  LEWIS  [Ref.  3,  ch.  ,6]  and 
LEWIS  [Ref.  9].  The  main  physical  reason  for  assuming 
independent  increments,  and  therefore  a  Poisson  process,  is 
that  the  operational  process  is  the  superposition  of  many 
individual  event  streams.  For  instance,  in  a  computer 
center  equipped  with  several  interactive  time-sharing 
terminals,  the  event  stream  of  users  at  each  specific 
teriinal  night  be  assumed  to  be  an  arbitrary  point  process 
with  a  certain  intensity  function.  The  event  stream  seen  by 
the  central  processing  unit  is  then  the  sum  total  (or 
superposition)   of   the   event   streams   of   the   individual 
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terminals  and,  if  there  are  enough  terminals,  should  be 
approximately  Poisson.  This  property  of  superposition  of 
intensity  functions  is  discussed  in  the  next  paragraph. 


A   Cal 1 s  per  hour 
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D.   DECOMPOSITION  AND  SUPERPOSITION  OF  POISSON  PROCESSES. 


To   obtain   a   Poisson   process   with  intensity  function 
A  (t)  =  X-,  (t)  +  \j  (t)  ,   we    may    superpose   two    Pcisson 
processes,  each  of  intensity  A-,(t)  and  /\2(t)» 

We  may  decompose  the  intensity  function  of  any  event 
stream  into  two  or  more  component  event  streams.  However  if 
we  then  superpose  the  component  streams,  we  will  recover  the 
original  type  process  only  if  we  started  with  a  Pcisson 
process.  For  example  begin  with  a  renewal  process  with 
intensity  v  (t)  =  and  let  v  =v,  +  v_  .  If  two  renewal 
processes  with  intensities  v,  and  v_  are  superposed,  the 
resulting  process  is  not  a  renewal  process. 

This  ucigue  property  may  be  exploited  when  simulating 
Poisson  processes.  It  permits  the  partitioning  cf  the 
intensity  function  to  take  advantage  of  any  special 
properties  cf  its  component  parts.  This  is  the  basis  for 
the  methcd  used  in  the  Poisson-Decomposition  and  Gap 
Statistic  Algorithm  discussed  in  later  sections. 


E.   TWO  EASIC  METHODS  CF  GENERATING  A  NON- HOMOGENEO OS 
PCISSON  EBOCESS 


1 •   lilSzJcale  Transformation  Algorithm 


Consider  the  non-homogeneous  Poisson  process  with 
intensity  function  (t)  >  0,  0<t<t0,  on  the  interval 
(0,t-.].    The   integral   of   the   intensity  function  is  then 
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A  (t)  and  the  number  of  events  in  the  interval  is  a   Pcisson 
random   variable   with    parameter  A(tQ)  -  A  (0)  =  y    (or 
A(tQ)  =  yQ  since   A  (0)  =  /°  X(t)  dt  =  0).   Now  if  T* ,  . . . ,  T* 
are   events   in   a   unit  homogeneous  Poisson  process  (i.e.  a 
Poisson    process    with    constant    intensity    function 

X  •  (t)  =  X1  =  1)  then  A"1  (T*)  , 
non-homogenecus  Poisson  process 


X  '  (t)  =  X1  =  1)  then  A-1  (T*)  ,  ...,A  1  (T*)  are  events  in  the 

1  n 


This  result  gives  a  procedure  for  simulating  a 
ncn-homogeneous  Poisson  process,  starting  with  a  homogeneous 
Poisson  process,  which  is  analogous  to  the  probability 
integral  transform  method  of  producing  random  samples  from  a 
continuous  distribution  with  distribution  function  F(x)  when 
starting  with  uniformly  distributed  random  variables.  The 
latter  method  is  essentially  that  if  the  inverse  of  F ( x)  can 
be  found,  then  by   generating   uniform   (0,1)   variates   u  , 

. .  . ,  u   the  values 

n 

X,  =  F~  (u,  )  ,  ...  ,  X   =  F~  (u  ) 

11  n        n 

comprise  a  sample  of  variates  from  the  desired  distribution. 

The  right-continuous  monotone  increasing  function 
A  (t)  describing  the  ncn-homogeneous  Pcisson  process  on  the 
interval  (0,tQ]  can  be  thought  of  as  a  distribution  function 
which   has   been   "scaled"   by   the    factor  u q  .     (Since 

A  (t Q)  =  UQ  '   then    A  ("t  J /y  '    1#   thus   A<t)/y0  is  a  v*lid 
distribution  function  on  (0,t  ].) 

Tc  implement  the  time-scale  transformation  procedure 
one  can  use  the  following  basic  result  for  homogeneous 
Pcisson  processes:  given  that  n  events  have  occurred  in  a 
homogeneous  Pcisson  process  over  a  fixed  interval  (C,tQ], 
those  events  are  uniformly  distributed  on  the  interval.  A 
prcof  of  this  property  is  given  in  PAHZSN  [Ref.  14,  p.  140], 
Therefore,  if  events  are  generated  as  a  unit  Poisson  process 
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on  an  interval  of  length  Ug  by  first  obtaining  a 
Poisson  (uQ)  random  variate  n,  and  then  letting  the  order 
statistics  from  a  random  sample  of  uniform  (0 ry q )  variates 
be  the  times  to  events  in  the  unit  Poisson  process,  and  the 
times  to  these  events  are  then  transformed  by  the  inverse  of 
the  integrated  intensity  function,  i.e.  A  (•) ,  the  effect 
is  the  same  as  that  obtained  from  the  probability  integral 
transform  method.  Figure  3  illustrates  this  method 
graphically  and  also  provides  a  flow  chart.  Note  the 
difference  however  between  this  procedure  and  the 
probability  integral  transform  procedure.  In  generating  a 
unit  Poisson  process  by  this  method  we  need  a  sample  whose 
size  is  random  (and  could  be  zero)  i.e.  Poisson  (  un)  •  The 
probability  integral  transform  method  simply  transforms  a 
fixed  number  cf  uniform  (0,1)  variates  into  variates  from 
some  other  distribution. 

In  Appendix  A  it  is  proved  that  scaling  both  the 
distribution  function  F (x)  and  the  interval  (0,1)  ty  the 
same  factcr  dees  not  affect  the  validity  of  the  probability 
integral  transform  method. 

The  unit  homogeneous  Poisson  process  may  also  be 
generated  by  adding  a  sequence  of  unit  exponential  variates 
(variates  frcm  an  exponentially  distributed  random  variable 
with  mean  =  1)  until  the  sequence  of  partial  sums  of  the 
random  variables  first  exceeds  u«  (see  Appendix  D) .  This 
accomplishes  two  things.  First  it  provides  an  ordered 
sequence  cf  events  from  a  homogeneous  Poisson  process  on  the 
interval  (0,Uq].  Second,  it  determines  a  realization  of  the 
Poisson  random  variable  Nt  with  parameter  A(tg)  =  Pq  .  Note 
that  given  n,  the  times  to  events  are  uniform  order 
statistics. 
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The  difference  between  these  methods  for  generating 
a  homogeneous  Poisson  process  should  also  be  noted.  The 
first  method  requires  a  Poisson  variate  and  ordering  of  that 
number  of  uniform  random  variables.  The  second  requires 
generation  of  independent  exponential  variates.  The  second 
method  is  probably  most  basic  in  that  it  requires  only 
exponentially  distributed  random  variables,  and  these  can  be 
obtained  from  uniform  variates  by  the  inverse  probability 
integral  transform,  which  is  just  a  logarithm.  The  method 
is,  however,  not  always  the  most  efficient. 

2 .    Conditioning   and   Order  Statistics  Algorithm  for  a 
Poisson  Process 


This  method  requires  the  result  of  a  theorem 
sketched  by  LEWIS  and  COX  [Bef.  3,  ch.  2]  and  restated  here 
for  convenience  and  continuity;  it  is  an  extension  of  the 
result  on  conditioning  in  a  homogeneous  Poisson  process 
which  results  in  a  conditional  uniform  distribution  of  the 
times  to  events. 

Theorem  J:   Assume  that  a  non-homogeneous  Poisson  process  is 

observed  for  a  fixed  time  (0,tg],   so   that   the   number   of 

events   in   (0,tg],   N.  ,   has   a   Poisson  distribution  with 

parameter   A(tn)  -   A(0)  =  u„    .  Then  if  T,,  .  .  .  ,T   denote 
0  0  1  n 

times-to-events  for  the  process  in  (0,t  ],  and  if  N,  =  n, 
conditional  on  having  observed  n  (>0)  events  in  (0 , t  ] ,  the 
T^'s  are  distributed  as  order  statistics  from  a  sample  with 
distribution  function 


FfH    ;\(t)  -  A(0) 
(    '         A(tQ)-  A(0) 


This    theorem    reduces    the    problem    of   simulating 
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non-homogeneous  Poisson  process   to  that   of  generating  a 

Pcisson  number  of  order  statistics  from  a  fixed  distribution 

function.   That  is,  given   an   intensity   function   over   an 

interval   (a,b],   whose  definite  integral  A (t)  =  /  X(s)ds  is 

a 

bounded   and   right-continuous   on      the      interval,      an      ordered 
sample,    ficm   a   population    with   distribution   function 


F(t) 


A(t)    "   A(a)  a    <    t        b  (1) 

A(b)    -   A(a) 


will  yield  the  desired  non-homogeneous  Poisson  process 
defined  by  the  intensity  function  X(t)  on  (a,b].  For 
simplicity  the  interval  will  hereafter  be  assumed  to  have 
its  left  end  point  at  zero  (a  =  0)  and  its  right  end  point 
at  some  arbitrary,  but  fixed  point  t  ,  (b  =  t  ) .  Using  this 
(0,t  ]  interval  results  in  no  loss  of  generality,  and  (1) 
becomes  identical  with  the  expression  in  the  theorem. 


Many  methods  exist  for  obtaining  the  necessary  order 
statistics.  The  inverse  integral  transform  explained  above, 
decomposition  of  the  density  function  (see  LEWIS  Ref .  8), 
or  the  rejection-acceptance  method  discussed  later  in 
Section  IV-D,  are  all  possibilities. 

For  the  family  of  intensity  functions  addressed  in 
this  thesis  the  non-homogeneous  Poisson  process  is  obtained 
by  a  combination  of  Poisson  decomposition,  an  algorithm  of 
LEWIS  and  SEFDLER  (Ref.  13]  for  obtaining  a  non-homogeneous 
Pcisson  process  with  log-linear  intensity  function,  and  the 
conditioning  and  order  statistics  theorem  given  above.  The 
procedure  involves  four  basic  steps. 

1.   The   intensity  function  is  decomposed  into  two 
components;  i.e.   A(t)  =  \(t)    +  A*  (t)  . 
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2.  A  gap  statistics  algorithm  is  used  to  generate 
a  ncn-homogeneous  Poisson  process  from  one  of  the 
components,  A  (t)  ,  which  is  chosen  such  that  it  has 
a  special  structure  (log-linear) . 

3.  A  rejection-acceptance  routine  is  used  to 
produce  a  sample  from  the  remaining  component, 
\*(t),  i.e.  a  Poisson  number  of  ordered  variates 
are  generated.  This  algorithm  is  described  in 
detail  in  Section  III-C.  This  sample,  when 
ordered,  becomes  a  Poisson  process  also. 

4.  The  Poisson  processes  of  the  component 
intensity  functions  are  then  superposed  to  produce 
the  desired  non-homogeneous  Poisson  process. 

Figures  4,  5  and-  6  (Section  III)  illustrate  the  four  general 
steps  of  this  procedure. 

Ncte:   Theorem  1  provides  a  second  way   to   generate 

the    unit   homogeneous   Poisson   process   required   by   the 

time-scale  transformation   algorithm   previously   described. 

First   we   may   generate   a   variate  from  the  Poisson  random 

variable   N+.   with   parameter  yn  ,    say    N.   =  n.     Then 
t0        r  0       *  t0 

conditional   upon   N.  =  n,  n  uniform  order  statistics  can  be 

?0 
generated   en   the   interval.    For   reasons   explained    in 

Appendix   D,   this   second   method   was  used  in  the  computer 

program   implementation   of   the   time-scale   transformation 

method. 
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F.       RATIONALE    FOR    SELECTION    OF    DEGREE-TWO    EXPONENTIAL 
POLYNOMIAL    INTENSITY    FUNCTION 


The    intensity    function    identifying      the      non-homogeneous 
Poisson    process   investigated    in   this   paper    is  of   the    form, 

2 

A(t)    =   exp(aQ    +    at   +   ct^t    ) 

where  a.,  a,  and  a9  are  real  constants. 

This  intensity  function  was  selected  for  three  reasons. 
First,  an  intensity  function  must  always  be  positive  (or 
zero)  if  it  is  to  be  meaningful.  The  above  function  is 
non-negative  for  all  values  of  a  /  a,  and  a  •  Secondly, 
this  intensity  function,  by  proper  choice  of  constants,  can 
be  used  tc  represent  the  four  different  types  of  event 
streams  mentioned  previously  in  Section  I.  And  finally,  the 
selection  of  this  intensity  function  leads  to  simple 
statistical  procedures;  (for  details,  see  LEWIS  Ref.  9, 
p. 30-34)  . 
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III.   COMPETING  ALGORITHMS 


A.   GENEEAL 


Assuming  an  intensity  function  of  the  form 


2 
A(t)  =  exp(an  +  a,  +  at  ) 

three  algorithms  for  generating  the  corresponding 
ncn-homogeneous  Poisson  process  are  discussed.  These 
algorithms  are  based  on  the  two  general  methods  presented  in 
Section  II-E,  and  the  decomposition  and  superposition 
property  of  Poisson  processes. 


B.   TIME-SCALE  TRANSFORMATION  ALGORITHM  (ALGORITHM  A) 


1  .   S t ejD  One 

Ey  definition  of  a  non-homogeneous  Poisson  process, 
the  total  number  of  events  observed  over  a  fixed  interval 
(0,tg)  is  itself  a  Poisson  distributed  random  variable,  Nt  , 

t  A  0 

with  parameter  Uq  =  Jq  A(t)dt.  The  first  step  cf  the 
algorithm  is  to  determine  the  value  of  the  parameter  Ug  . 
Although  an  explicit,  closed-form  expression  for  the  above 
integral  cannot,  be  found,  a  series  representation  does 
exist.     Except    for    a    constant   factor,   this   series 
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representation  assumes  the  form  of  the  error  function  or  of 
Dawson's  integral.  A  negative  value  for  the  coefficient  of 
the  second  degree  term,   a  $    yields  the  error  function  form: 


Uft  -  K(  2-  !    2   e"u   du  -  i-  /  1   e"u   du 


t„     2       „   t,       _„2 
-  K  [  -£-  / 

/?0  ^  0 


whereas      a      positive      value      for  results   in   the    Dawson's 

integral   fori: 

„„.«»(*;?    j^d.)     -   KtJ  (t-    I'   S^   du) 


In  the  above  expressions  K,  t  ,  and  t_  are  uniguely 
determined  by  the  coefficients  an,  a,,  a-  and  the  end  points 
of  the  interval  over  which  the  intensity  function  applies. 
A  detailed  derivation  of  above  relationships  is  given  in 
Appendix  C. 

Evaluation  of  the  error  function  is  a  FORTRAN 
supplied  procedure  and  requires  only  that  the  proper 
arguments  be  calculated  and  provided  to  the  FORTRAN  FUNCTION 
ERF  or  EERF  [Ref.  15].  Evaluation  of  Dawson's  integral  is 
best  accomplished  through  use  of  the  IMSL  (International 
Mathematical  and  Statistical  Libraries,  Inc.)  FUNCTION  MMDAW 
[Ref.  6].  The  accuracy  of  the  function  values  calculated  by 
these  routines  is  limited  only  by  the  precision 
characteristics  of  the  computer. 


2  .   S t e_p  Two 

Cnce  the  parameter  of  the   Poisson   random   variable 
N     is  determined,  a  realization  on  that  random  variable  is 


t 


0 
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required.  (The  approach  is  somewhat  backwards  since  it 
first  determines  hew  many  events  occurred  over  the  interval. 
It  then  distributes  that  fixed  number  of  events  over  the 
interval  in  accordance  with  the  non-homogeneous  Pcisson 
process  described  by  the  intensity  function.  The  importance 
of  Theorem  1  now  becomes  evident  since  it  assures  the 
validity  of  such  a  procedure.)  Generation  of  Poisson 
variates,  especially  those  with  large  parameter  values,  is  a 
complex  procedure  in  itself  if  efficiency  in  terms  of 
computer  time  and  memory  requirements  is  desired.  This 
problem  is  discussed  later  in  Section  IV-B.  For  the 
present,    assume    that   the   requisite   variate   has   been 


produced,  i.e.  N    =  n. 

fc0 

3 •   J t e£  Three 

Given  that  n  events  have  occurred  over  the  interval 
(0,tg]  we  then  distribute  n  events  along  an  interval  of 
length  Li  in  accordance  with  a  homogeneous  Poisson  process. 
Since  events  in  a  homogeneous  Pcisson  process  are  uniformly 
distributed  over  an  interval  (given  that  n  events  have 
occurred),  this  step  merely  requires  that  n  uniform  (0,1) 
variates  be  generated,  ordered. from  lowest  to  highest,  and 
then  each  multiplied  by  the  factor  u_.  The  values  in  this 
n-element  vector,  (u,'  ,  u '  ,  ...,  u')  ,  correspond  to  the 
pcints  plotted  on  the  vertical  axis  in  Figure  3. 

4 .   s t ep  Four 

Each   event   in  the  homogeneous  Poisson  process  must 
be  transformed   by   the   inverse   of   the   integral   of   the 

intensity   function.   Letting  /   A  (s) ds  =  A  (t) ,  the  inverse 

i  ° 

A"1  (•)  applied  to  each  event   in   the   homogeneous   Pcisson 

process,    will    produce    a   corresponding   event   in   the 
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non-homogeneous   Poisscn    process,    i.e.      A   (u,)  =  t,  , 

-1  • 

A   (u_)  =  t  ,   etc.    The   difficulty   is   that   since   the 

integral  of  this  specific  intensity  function  cannot   usually 

be   explicitly   expressed,   the   form  of  its  inverse  usually 

eludes  any  convenient  computational  formula  expression.   The 

unique   pcsiticn   on   (0,t-]  the  inverse  determines  for  each 

input  value  can  be  found  to  any  degree  of  accuracy   desired, 

by   iterative,  numerical  methods.   The  Newton-Raphson  method 

is   easily   employed   and   very   efficient   in   the   present 

scenario.    Its  implementation  is  explained  in  Section  IV-C. 

Since  the  function   A  (t)  is   strictly   monotone   increasing, 

the   inverse  function   A   (u)  applied  to  an  ordered  sequence 

of  input  values  results  in  an  ordered   sequence   of   output 

values.    Therefore,  t,  ,  t„,  ...,  t   are  the  times  cf  events 

12        n 

in  the  non-hcmogeneous  Poisson  process  and  the  algorithm   is 
complete . 


C.   TIME-SCALE  TRANSFORMATION  ALGORITHM,  ALTERNATE 
(ALGORITHM  A') 


An  alternative  approach  to  the  time-scale  transformation 
method  described  above  is  to  generate  the  required 
homogeneous  Poisson  process  by  using  the  fact  that  in  this 
process  the  random  times  between  events  are  independently 
exponentially  distributed.  Thus  one  generates  unit 
exponential  variates  until  their  sum  exceeds  n  .  The 
partial  sums  give  the  times  to  events  and  the  number  of 
partial  sums  less  than  or  equal  to  u  is  a  Poisson  (,,  ) 
variate.  Note  that  the  Poisson  variate  comes  out  as  a 
by-product  in  this  procedure  rather  than  as  a  pre-product  as 
in  Step  Two  of  Algorithm  A  above. 

Although  this  method  combines  Step  Two  and  Step  Three  of 
Algorithm  A  into  a  single  procedure,  it  is   not   necessarily 
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the  best  method.  Because  it  requires  the  use  of  the 
additive  method  of  Poisson  variate  generation,  it  becomes 
inefficient  for  Poisson  processes  with  many  events  (see 
Appendix  D) .  For  this  reason  Algorithm  A  instead  of 
Algorithm  A'  was  used  when  implementing  the  time-scale 
transformation  method  into  a  computer  program. 


D.   POISSCN-DECOMPOSITION  AND  GAP  STATISTIC  ALGORITHM 
(ALGORITHM  B) 


It   is  recommended  that  the  reader  refer  to  Figures  4,  5 
and  6  for  a  better  understanding  of  the  following  steps.) 

1 .   S t e£  One 

The  Poisson- decomposit ion  and  gap  statistic 
algorithm  begins  with  an  examination  of  the  coefficients  of 
the  intensity  function.  By  doing  so,  the  intensity  function 
is  categorized  into  one  of  six  possible  configurations. 
These  six  cases  are  discussed  in  LEWIS  and  SHEDLEH 
[Ref.  11].  Examples  of  each  case  are  illustrated  in  Figures 
7  through  12  located  at  the  end  of  this  section. 

2  .   St e£  Two 

a.  If  the  intensity  function  \  (t)  is  monotone 
increasing  or  monotone  descreasing  over  the  interval  (Cases 
I,  II,  IV  and  V;  see  Figures  7,  8,  10  and  11)  the  intensity 
function  is  decomposed  into  two  separate  intensity  functions 
over  the  same  interval.  The  resulting  intensity  functions 
are  of  the  fcrm; 

A(t)  =  exp(YQ  +  Yjt) 
and 
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A*(t)  =  X(t)  -  A(t)  =  exp(aQ  +  a1t  +  a2t2)  -  exp(yQ  +  y^t) 
It  is  clear  that   A (t)  +  A*(t)  =  X(t). 


b.  If  either  of  the  two  cases  not  covered  by  2a 
abcve  occurs,  A(t)  will  be  monotone  increasing  (decreasing) 
on  the  subinterval  (0 , b  ]  and  monotone  decreasing 
(increasing)  on  the  complementary  subinterval  (bftg]  (see 
Figures  3  and  6) ,  where  b  is  a  unique  point  within  the 
interval  at  which  X (t)  has  a  maximum  (minimum)  value.  By 
dividing  the  interval  properly  into  two  disjoint,  contiguous 
subintervals,  each  subinterval  may  be  treated  as  explained 
in  2a.  Subsequent  steps  are  applied  to  each  of  the  two 
intervals  separately  and  the  results  combined. 

3 •   St e£  Three 

An  efficient  algorithm  for  generating  a 
ncn-homogeneous  Poisscn  process  with  a  log-linear  intensity 
function  (i.e.  A ( t)  =  exp  ($  ♦  3,t))  is  presented  by  LEWIS 
and  SHZELEE  [Hef.  13].  This  algorithm  generates  the 
ncn-homogenecus  Poisson  process  through  the  use  of  gap 
statistics.  (A  comparison  of  the  gap  statistic  technique 
with  the  conventional  integral  transform  technique  is 
discussed  in  Appendix  C.)  By  judicious  selection  cf  the 
coefficients  cf  the  log- linear  intensity  function,  most  of 
the  total  area  under  the  original  intensity  function  A  (t) 
will  be  contained  under  the  function  A (t) .  Therefore  most 
of  the  events  in  the  non- homogeneous  Poisson  process  with 
intensity  function  A  (t)  can  be  accounted  for  by  employing 
the  gap  statistics  algorithm  on  the  intensity  function 
X  (t)  . 
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Step  1  -  Categorize  intensity  function 


Step  2 


1  2 

Decompose     \(t)     into     A(t)    (log-linear)  and 
A*(t)   =  A(t)   -  A(t). 


A 


6" 


4-. 


A*(t) 


Ficure    U    - 


DIAGRAM    OF    POISSCN-DECOMPOSITION     AXID    GAP 
STATISTIC    ALGORITHM 
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ti  t2  HH       hH     t7  t 


8 


i  i  r 

0  1  2 

Step  3  -  Generate  events  T.  from  log-linear  intensity 

function,  A_(t),  using  gap  statistic  algorithm. 
Events  will  be  ordered. 


s2       » 


Step  4  -  Generate  events  s.  from  intensity  function 

X*(t)  =  x(t)  -  _X(t)   using  rejection  algorithm, 
Events  will  not  be  ordered. 


S(D      S(2) 


0  1 

Step  5a  -  Order  events  from  Step  4 


Tl   T2     T3T4  T5  T6         T7  T8  T9  T10 

1111  111  II  I 


0  1  2 

Step  5b  -  Merge   (superpose)  events   from  Steps  3  and  5a 
Result   is  event  stream     T, ,   T   ,    ..,      from 
\(t)   =   X(t)   +  A*(t). 

Figure    5    -       DIAGRAM     (CONTINUED) 
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Streams 


Figure    6    - 


FLOW    DIAGRAM    OF    POISSON-DECOMPOSITIO M    AND    GAP 
STATISTIC    ALGORITHM 
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4  .   S tejj  Four 


All  that  remains  to  be  done  is  to  generate  an 
ordered  sample  of  some  given  size,  on  the  interval  (0,tn], 
frcm  the  remaining  component  of  the  original  intensity 
function,  i.e.   X*  (t)  .. 


The  number  of  events  in  the  interval   is   a   Poisson 

random   variable   N!     with   parameter     u '  =   /  0\*(t)dt. 

z0  0 

Once  a  realization  on  this  random  variable  is  obtained,  i.e. 

N*   =  n'f   Theorem   1   may  be  invoked.   Since    A*(t)/ui  is 
t0  0 

more  easily  evaluated  than  its  indefinite  integral,  the 
rejection  technique  (explained  in  Section  I7-D)  is  used  to 
generate  the  n*  required  variates.  The  rejection  technique 
is  not,  in  general,  always  an  efficient  method  for  variate 
generation  unless  great  care  is  taken.  Yet  in  this 
deccmpositicn  scenario,  it  will  be  used  to  generate  only  a 
small  percent  of  the  total  events  required  by  the  original 
intensity  function  A(t).  The  majority  of  the  events  will 
be  generated  by  the  efficient  gap  statistics  algorithm.  The 
efficiency  gains  should  more  than  compensate  for  any 
efficiency  lcsses  due  to  the  use  of  the  rejection  technique. 

5.   Step  Five 


The  events  produced  by  Step  3  will  be  in  order  on 
the  interval  (0,tg].  The  events  produced  in  Step  4  will  not 
be  in  order  en  (0,tg].  By  ordering  the  events  from  Step  4, 
(which  are  few  in  number  compared  to  the  total  number  of 
events  in  the  non-homogeneous  Poisson  process)  it  is 
possible  tc  superpose  the  two  ordered  event  streams.  The 
merged  event  streams  produce  a  new  event  stream  from  the 
original  intensity  function    X  (t) . 
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figure    7    -       SAMPLE    INTENSITY    FUNCTION    -    CASE    I 
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Figure  9  -   SAMPLE  INTENSITY  FUNCTION  -  CAS3  III 
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IV.        ALGORITHM    IMPLEMENTATION 


A.        GENERAL 


This  section  explains  the  several  specific  techniques 
that  were  applied  to  implement  the  two  competing  algorithms 
(Algorithm  A  and  Algorithm  B)  into  FORTRAN  computer 
programs.  Detailed  discussion  of  the  various  subprograms  is 
avoided  since  References  11  and  13  and  attached  program 
listings  provide  such  information.  The  Algorithms  A  and  B 
have  some  subprogram  requirements  in  common  while  other 
subprograms  are  unique  to  one  algorithm  or  the  other.  This 
section  will  discuss  first  those  requirements  common  to  both 
algorithms,  then  those  needed  only  by  the  time-scale 
transformation  algorithm  (Algorithm  A)  and  finally  those 
unique  to  the  Poisson-decompositio n  and  gap  statistic 
algorithm  (Algorithm  B) .  Hereafter  differentiation  between 
the  algorithms  and  the  FORTRAN  computer  program 
implementation  of  the  algorithms  will  not  always  be  made. 
The  meaning  will  be  clear  from  the  context. 


B.   COMMON  REQUIREMENTS 


1 ♦   Integration  of  a  Degree-Two   Exponential   P  cljromial 
Function  over  a  Fixed  Interval 


Algorithm  A  requires   that   the   intensity   function 
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X  (t)  be  integrated  ever  a  fixed  interval  in  order  to 
determine  the  value  of  the  parameter  of  the  Poisson  random 
variable  that  governs  the  number  of  events  that  are  observed 
to  occur  in  the  non-homogeneous  Poisson  process.  Algorithm  B 
reguires  that  the  intensity  function  X*(t)  =  X(t)  -  X  (t)  , 
be  integrated  over  a  fixed  interval  for  the  same  reason. 
Since 

b  b  b 

/   X*(t)dt  =  /   X(t)dt  -  /   X(t)dt 
a  a  a 


and   X(t)  has  an   explicit   expression   for   its   indefinite 
integral ,  i. e. 


/   X(t)dt  =  exp(yn) 


exp(Y-jt) 


the  problem  fcr  both  algorithms  is  reduced  to  computing  the 
value  for 

b  b  2 

/   X(t)dt  =  /   exp(aQ  +  a..t  +  c^t  )  dt  . 
a  a 

S0ERO0TINE  HELP  employs  IMSL  FUNCTION  MMDAW  or  the  FORTRAN 
supplied  procedure  DERFr  as  appropriate,  to  perform  this 
calculation.  Section  III-3  and  Appendix  B  discuss  how  this 
computation  could  also  be  made  using  a  convergent  series 
representation. 

2 •    Ge neration   of  a   Poisson   Varia te   with   a   Given 
Parameter 


Both  candidate  algorithms  require  at  least  one 
realization  on  a  Poisson  distributed  random  variable  with  a 
given   parameter.    The  value  assumed  by  the  random  variable 
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is  the  number  of  events  observed  in  a  specific  interval  of 
time  over  which  the  occurrence  of  events  is  governed  by  a 
given  intensity  function  (i.e.  X  (t)  for  Algorithm  A  and 
X  *(t)  fcr  Algorithm  B)  .  The  nature  of  most  real  event 
streams  which  lend  themselves  to  analysis  using  a 
non-homogeneous  Poisson  process  model  is  that  they  consist 
of  a  large  tctal  number  of  events.  This  will  be  the  case 
either  if  a  dense  process  is  observed  over  an  interval  of 
shcrt  or  moderate  length  or  if  a  "sparse"  process  is 
observed  over  an  extended  interval  (see  the  arrivals  at  an 
intensive  care  unit  given  in  LEWIS  [  Ref .  9]).  The  point  is 
that  both  algorithms  must  be  flexible  enough  to  generate 
non-homogeneous  Poisscn  processes  that  result  in  high 
numbers  cf  events  occurring.  Since  the  number  of  events 
occurring  over  any  interval  in  such  a  process  is  a  Poisson 
random  variable,  it  becomes  necessary  to  be  able  to  generate 
a  variate  from  a  Poisson  distribution  with  a  large 
parameter,  i.e.  large  mean.  The  two  most  theoretically 
straightforward  algorithms  for  generating  Poisson 
distributed  variates  (the  additive  and  multiplicative 
methods)  become  computationally  cumbersome  as  the  parameter 
of  the  random  variable  increases.  Both  the  additive  and 
multiplicative  algorithms  and  the  practical  deficiencies  of 
each  are  discussed  in  Appendix  D. 

The  technique  selected  to  deliver  a  Poisson  variate 
on  demand  to  both  Algorithm  A  and  Algorithem  B  is  the  Gamma 
method.  It  is  developed  and  explained  in  an  unpublished 
book  by  AHEENS  and  DIETER  [Ref.  1].  A  paraphrased  account 
of  the  development  of  this  algorithm  is  included  in 
Appendix  E.  The  main  advantage  of  AHRENS  and  DIETER'S  Gamma 
method  is  that  whereas  the  computer  time  required  for  the 
additive  and  multiplicative  methods  is  proportional  to  the 
parameter  of  the  Poisson  distribution  being  sampled,  the 
computer  time  required  by  the  Gamma  method  is  proportional 
to  the  logarithm   of   the   parameter.   The   SUBROUTINE   PVAR 
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employs  the  Gamma  algorithm  to  return  a  Poisson  variate  when 
given  a  parameter  as  an  input. 

3«   IX£Ht  Storage 

Any  algorithm  which  simulates  a  non-homogeneous 
Poisson  process  must  have  some  mechanism  to  provide  the  user 
with  information  that  completely  describes  the 
realization (s)  of  the  process  simulated.  Since  the  specific 
arrangement  of  the  stream  of  events  over  the  interval  is  the 
information  cf  interest  to  the  analyst,  the  location,  or 
time  of  occurrence,  of  each  single  event  on  the  interval 
must  be  stored.  Eguivalently ,  the  spacings  between  events 
would  completely  describe  the  realizations  on  a 
non-homogeneous  Poisson  process.  Thus  event  spacing 
information  rather  than  event  location  information  could  be 
stored.  (Programs  written  for  candidate  algorithms  A  and  B 
both  provide  the  option  for  the  user  to  demand  either  event 
location  or  event  spacing  information.)  When  using  either 
algorithm,  an  array  large  enough  to  hold  location  or  spacing 
data  for  each  event  generated  by  the  algorithm  must  be 
created.  Since  the  number  of  events  observed  in  any 
realization  of  a  non-homogeneous  Poisson  process  is  itself  a 
random  variable,  the  precise  size  of  the  array  cannot  be 
determined  a  priori.  If  the  programs  are  to  have  any  value 
for  general  application,  they  must  be  able  to  accept 
intensity  function  parameters  which  will  demand  large 
numbers  of  events  when  simulated.  Using  the  somewhat 
arbitrary  assumption  that  an  event  stream  with  an  average  of 
4500  events  is  sufficient  for  most  simulation  scenarios,  a 
fixed  array  with  a  capacity  for  5000  events  is  used  in  the 
programs  implementing  both  of  the  algorithms.  If  the  value 
of  the  intensity  function  integrated  over  the  interval  is  as 
high  as  the  maximum  limit  of  4500  (i.e.  the  number  of  events 
in   the   interval   is   Poisson  distributed  with  mean  =  4500) 
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then  the  array  of  length  5000  allows  the  Poisson  random 
variable  to  exceed  its  mean  by  7.45  standard  deviations 
before  an  array  overflow  is  encountered.  This  highly 
unlikely  event  is  of  such  rare  occurrence  (less  than  one 
chance  in  a  billion)  that  it  may  be  disregarded.  However, 
should  it  occur,  the  programs  will  reinitialize  themselves 
and  generate  a  new  Poisson  process.  Also  an  error  indicator 
is  incremented.  This  error  indicator  (IER)  may  be  written 
on  demand.  Its  value  is  the  number  of  times  the  program  was 
forced  tc  abort  generating  a  Poisson  process,  reinitialize 
itself  and  start  again. 

A  50C0  element  array  is  small  enough  to  avoid  its 
being  an  undue  memory  requirement  burden  on  most  operating 
systems,  yet  large  enough  to  accommodate  most 
non-homogenecus  Poisson  processes  of  interest.  Its  choice, 
though  arbitrary,  was  based  upon  the  above  two 
considerations. 


C.   SPECIAL  REQUIREMENTS  OF  THE  TIME-SCALE  TRANSFORMATION 
ALGORITHM  (ALGORITHM  A) 


1  •   Dnif orm  Variates 

As  explained  in  Section  III-B,  once  the  number  of 
events  observed  in  the  non-homogeneous  Poisson  process  is 
established  (i.e.  N  =  n)  ,  it  is  necessary  to  construct  a 
homogeneous  Poisson  process  consisting  of  n  events  over  an 
interval  of  length  yQ  units.  This  is  easily  done  by 
generating  n  uniformly  distributed  variates  on  the  interval 
(0,1)  and  then  scaling  each  by  the  factor  yn  .  The  LLRANDOM 
computer  library  package  developed  by  LEWIS  and  LEARMONTH 
[Ref.  12]  is  a  very  efficient  source  of  such  variates.   Once 
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an  array  of  n  uniform  (0,1)  variates  is  obtained  from 
LLRANDOM,  each  element  of  the  array  must  be  multiplied  by 
th€  appropriate  scaling  factor.  It  is  these  resulting 
numbers  which  must  be  acted  upon  by  the  inverse  of  the 
integrated  intensity  function  to  yield  the  location  of 
events  in  the  non-homogeneous  Poisson  process  on  the 
original  interval  (0,tQ]. 

2  •   Sorting  of  Events 

To  be  of  much  practical  value  a  simulation  routine 
for  a  ncn-hcmogeneous  Poisson  process  should  order  the 
events  in  the  interval  from  first  to  last.  In  Algorithm  A, 
this  ordering  may  be  done  before  applying  the  inverse 
integrated  intensity  function  to  the  elements  in  the 
homogeneous  Poisson  process.  The  monotonic  nature  of  the 
integrated  intensity  function  maintains  the  relative  order 
of  all  elements  after  they  have  been  transformed  by  the 
inverse  function  (see  Figure  3) .  Of  course  the  ordering 
could  be  dene  after  the  transformations  also.  In 
implementing  Algorithm  A,  the  uniform  variates  were  scaled 
by  the  factor  UfJ  then  ordered,  from  lowest  to  highest, 
before  being  transformed  by  the  inverse  of  the  integrated 
intensity  function. 

Ordering  of  large  arrays  of  numbers  is  a  time 
consuming  operation  on  the  computer.  There  are  many 
ordering  algorithms  of  varying  degrees  of  sophistication  and 
efficiency.  Because  ordering  is  unavoidable  when  using 
Algorithm  A,  selection  of  an  efficient  ordering  routine  is 
important.  The  ordering  algorithm  used  in  this 
implementation  was  the  W.  R.  CHURCH  computer  center  library 
routine  FXSOBT  which  employs  Singleton's  version  of  the 
partition  exchange  sort.  (A  program  listing  of  PXSCRT  is 
provided  in   the   computer   center's   catalogue   of   library 
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routines.)  The  PXSORT  routine  appears  to  be  the  most 
efficient  ordering  routine  readily  available  for  the 
purpose.  (It  is  acknowledged  that  more  efficient  routines 
specifically  tailored  to  the  problem  of  ordering  uniform 
random  variates  may  possibly  have  been  developed  that  would 
improve  the  overall  efficiency  of  the  program  implementing 
Algorithm  A.) 

3.   Computation  of  the  Transformed  Values 

The  essence  of  Algorithm  A  is  the  application  of  the 
inverse  cf  the  integrated  intensity  function  to  each  event 
in  the  homogeneous  Poisson  process  on  the  interval  (0,uq]« 
As  mentioned  before,  the  intensity  function,  A(t) r  under 
consideration  does  not  yield  an  explicit  expression  for  the 
integrated  intensity  function,  A(t)  =  F(t) ,  that  must  be 
inverted.  (Note:  F (t)  is  a  "scaled"  distribution  function). 
Hence  neither  does  a  computationally  convenient  expression 
for  this  inverse  function  exist.  Numerical  methods  must  be 
employed  to  transform  the  position  of  each  event  on  the 
interval  (0,Uq]  in  the  homogeneous  process  to  its 
corresponding  position  on  the  interval  (0,tg]  over  which  the 
simulated  non-homogeneous  Poisson  process  is  to  be  produced. 

The  Newton-Raphson  technique  was  used  to  accomplish 
the  transformation.  This  technique  allows  for  iterative 
approximations  which  converge  to  the  true  transformed 
values.  (i.e.  t1#  t2,  ...,  tn;  where  t^  =  F  (u-))  The 
iterations  continue  for  each  value  u^  until  an  estimate  t' 
for    its    corresponding    t^    is    found   that   satisfies 

jF(t')  -  u.j  \    <    z    f    where  £  is   a   predetermined   tolerance 

-5 
level.   (The  specific  value  for   £  used  was  u«  x  10    .) 

In  the  Newton-Raphson  technique,  given  a  function 
h  (x)  ,   the  objective  is  to  find  the  solution.,  x*,  satisfying 
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the  equation  h  (x*)  =  0.   Letting  x-,  be  the  initial  estimate 

for   x*   (based   upon   some   reasonable   criteria)  the  next 

estimate  for  x*,  that  is,  x  2r  can   be   calculated  from   the 
fundamental  iteration  relationship, 


h(x  ) 
X,..,  =  x,  - 


k+1    ~k    h' (xk)   ' 

This  assumes  that  the  function  h (•)  and  its  derivative  h*(») 
can  be  evaluated  at  x^.  The  process  is  repeated  k  times 
until  |  h  (X]<)  |  <£  ;  or,  until  J  xk  -  xk+1 1  <  £  .  This 
procedure  is  nothing  more  than  using  the  first  two  terms  of 
the  Taylcr  series  expansion  of  the  function  h  (•)  to  locate 
the  x-intercept  of  the  line  tangent  to  h (x)  at  the  point 
h(xk).  That  intercept  value  becomes  the  next  approximation 
for  the  roct  of  h (x)  and  the  procedure  is  repeated.  A 
graphical  explanation  of  the  Newton-Raphson  method  is 
presented  in  Appendix  E. 

Applying  the  Newton-Raphson  method   to   the   present 

problem  requires  some  special  modifications.   The  problem  is 

to  find  a  value  t.  such  that  F  (t. )  =  u.,  where  u.  is   known, 

1  v  1'     i'        1 

i.e.   to   find   t.  =  F   (u.) .    Direct   apDlication   of   the 

inverse  F   (•)   is   impossible   since   its   functional   form 

defies   all   but   the  most  abstract  mathematical  expression. 

However  if  a  new  function  G.  (t)  =  F(t)  -  u.  is   defined   and 

if   its  root  t*,  determined  (i.e.  a  t*  such  that  G-  (t*)  =  0) 

then  F(t*)  -  u.  =0,  and  F(t*)  =  u.  .   Thus   t.  =  t*   is   the 
\   /     1     i  v   '     i  i 

value  desired. 

In  order  to  apply  the  Newton-Raphson  method  tc  the 
function  G  .  (t)  it  must  be  possible  to  calculate  both  G.  (t) 
and  its  derivative.  Since  G.  (t)  differs  from  F (t)  only  by  a 
constant,  the  problem  reduces  to  that  of  evaluating  F  (t) . 
This  has  already  been  done,  as  previously  explained  in 
Section    III-B,   and   merely   requires   the   assistance   of 
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SUEROUTINE   HELP.    Now   G '.  (t)  =  F1  (t)  ,    so    taking  the 

derivative   cf  G  .  (t)   returns   the  intensity  function  A  (t) 

which   is   easily   evaluated   for  any   t.     Clearly,  the 
Newton-Raphson  method  may  be  used. 

It  should  be  noted  here  that  for  each  u  .  there  is  a 
corresponding  function  G  . (t)  en  which  the  Newton-Raphson 
technigue  must  be  used.  Since  each  use  of  the 
Newton-Raphson  method  may  require  several  iterations  to 
arrive  at  a  value  t*  which  satisfies  the  tolerance 
criterion,  it  is  readily  evident  that  for  large  n, 
considerable  computational  effort  is  required  to  obtain  all 
the  values  t.  (i  =  1 ,  ...,  n) .  The  number  of  iterations 
required  to  arrive  at  a  suitable  approximation  for  each  t. 
is  highly  sensitive  to  the  accuracy  of  the  initial 
approximation  (designated  t  )  If  the  initial  approximation 
is  close  to  the  actual  t.,  the  Newton-Raphson  method  will 
converge  very  quickly.  If  the  initial  approximation  is 
poor,  convergence  could  be  much  slower.  The  procedure  used 
to  select  initial  approximations  for  each  t.  will  therefore 
have  a  profound  effect  upon  the  overall  efficiency  of 
Algorithm  A. 

Selection  cf  these  initial  approximations  is  dene  by 
partitioning  the  interval  (0,t  ]  into  equal  length  segments. 
The  number  of  segments  is  equal  to  min[10,  n/4].  The 
function  F  (•)  is  evaluated  at  the  end  points  of  each  segment 
and  these  end  points,  with  their  corresponding  function 
values  are  stored  in  an  array.  This  procedure  is  performed 
by  SUBROUTINE  BNCHMK.  See  Appendix  E  for  a  graphical 
representation. 

To  find  an  initial  approximation  for  the  t^  which 
corresponds  to  any  given  u.r  the  array  of  function  values  is 
searched  until  two  adjacent  function  values  which  bracket  u. 
are   located.    These   adjacent   function   values    uniquely 
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identify  the  segment  of  the  interval  (0,t  ]  in  which  the 
elusive  t.  is  located.  Either  end  point  cf  the  segment 
would  serve  as  a  good  initial  approximation  tj_  ^for  the 
Newton-Raphson  method.  However,  for  the  purpose  of  this 
implementation,  the  end  point  which  yielded  the  function 
value  closer  to  uj_  is  used  for  the  initial  approximation. 

The  decision  to  divide  the  interval  (0,tg]  into 
min  [10,  n/U  ]  segments  was  based  upon  empirical  results.  Of 
several  proposed  segmenting  schemes  that  one  which  resulted 
in  the  fewest  total  number  of  distribution  function 
evaluations  over  several  intensity  function/interval 
scenarios  was  chosen.  Higher  values  than  n/4  for  sparse 
event  streams  may  produce  better  results.  But  n/4  gave  good 
results  for  dense  streams  and  adequate  results  for  sparse 
streams.  It  appeared  to  be  the  best  compromise  as  a 
candidate  for  general  usage.  (The  number  of  function 
evaluations  of  Gj_(t)  for  each  t-  averaged  between  2.2  and 
2.7  for  this  partitioning  scheme.) 

Alternative  methods  would  be  to  use  for  the  initial 
approximation  for  t   one  of  the  following  values: 

*i-i  +  t/n 

Neither  cf  these  other  methods  were  attempted  so  it  is 
uncertain  whether  they  would  be  more  or  less  efficient  than 
the  method  which  was  used.  Efficiency  differences  among  the 
three  methods  are  probably  not  substantial  since  each  will 
usually  yield  a  good  first  approximation. 
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D.   SPECIAL  REQUIREMENTS  OF  THE  DECOMPOSITION  ALGORITHM 
(ALGORITHM  B) 


1 •   Intensity  Function  Categorization 

The  decomposition  routine  selected  depends  on  which 
of  six  possible  shape  categories  the  intensity  function 
A  (t)  falls  into  (c.f.  Figures  7  through  12)  .  An 
examination  cf  the  constants  C6q  ,  a^  and  o^  iQ  tne  intensity 
function  and  the  interval  (0,t0]  over  which  the 
non-homogeneous  Poisson  process  is  to  occur  will  uniquely 
identify  the  category  of  the  intensity  function.  A  thorough 
discussion  of  this  procedure  is  presented  in  Ref.  11,  and  is 
net  reproduced  here.  Implementation  of  the  category 
identification  procedure  requires  a  lengthy  sequence  of 
decision  statements  within  the  computer  program. 

2 •    Selection   of   the   Imbedded   Log-Linear   Intensity 
Function  (s) 

Once  the  intensity  function  has  been  categorized  it 
must  be  decomposed  in  accordance  with  a  complex  scheme 
described  in  Ref.  11.  The  objective  of  the  decomposition 
scheme  is  to  fit  a  leg-linear  curve  (or  curves)  completely 
underneath  A  (t)  on  the  interval  (i.e.  A  (t) <  A  (t)  ,  0<t<tQ) 
in  such  a  way  as  to  maximize  the  area  under  the  log-linear 
curve (s)  .  This  is  done  by  partitioning  of  the  interval 
(0,tQ]  into  subintervals  (0,b]  and  (brtQ]  if  necessary,  and 
then  by  proper  selection  of  the  coefficients  Yq  and  Yi  -or 
the  log-linear  f unction (s).  These  coefficients  are 
functions  cf  the  coefficients  a  ,  a   and  a   in  the   original 
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intensity  function  A  (t)  and  of  the  interval  (0,tQ].  As  the 
program  advances  through  the  categorizing  decision 
statements,  the  proper  coefficients  are  computed  fcr  the 
appropriate  leg-linear  function (s). 

3  •   Gap.  Statistics  Algorithm 

The  precursor  to  LEWIS  and  SHEDLER  [Ref.  11]  was  a 
paper  by  the  same  authors  [ Ref .  13]  proposing  a  gap 
statistics  algorithm  for  simulating  a  non-homogeneous 
Poisson    process    with   a   log-linear   intensity   function 

A  (s)  =  exp  (3q  +  BtS)  .  (Reference  11  reviews  this  algorithm 
in  detail.)  As  previously  noted,  Algorithm  B  divides  the 
intensity  function  X (t)  into  a  sum  of  two  intensity 
functions,   _A(t)   and    A*(t).    The  new  intensity  function 

A_(t)  is  chosen  to  take  on  the  form  of  a  log-linear  function 
so  that  the  cap  statistics  algorithm  may  be  used  to  generate 
a  stream  of  events  from  this  portion  of  the  original 
intensity  function  A (t)  .  The  SUBROUTINE  NHPP2  implements 
the  LEWIS  and  SHEDLER  gap  statistics  algorithm  for  the  input 
intensity  function  A_(t)  and  returns  an  appropriate  event 
stream  to  the  calling  program.  Since  the  use  of  the  gap 
statistics  algorithm  within  Algorithm  B  is  the  rationale  for 
suspecting  it  to  be  a  more  efficient  algorithm  than  the 
time-scale  transformation,  it  is  instructive  to  examine  the 
efficiencies  gained  in  using  the  gap  statistics  algorithm 
vice  a  time-scale  transformation  algorithm  on  a  leg-linear 
intensity  function.  This  question  is  addressed  in  Appendix 
C.  It  was  found  that  use  of  the  gap  statistics  algorithm 
resulted  in  a  reduction  in  computer  time  of  approximately 
50%  from  that  required  by  the  time-scale  transformation 
method. 
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1  •  2h§.    Beject ion  Routine 

The  gap  statistics  algorithm  has  efficiently 
produced  an  event  stream  from  a  non-homogeneous  Poisson 
process  defined  by  the  log-linear  intensity  function 
X  (t)  =  exp  (yn  ♦  y.  t)  .  But  the  process  to  be  simulated  has 
intensity  function  \  (t)  =  exp  (aQ  +  a,t  +  a  t2).  It  is 
therefore  necessary  to  superpose  an  event  stream  frcm  the 
intensity  function 

X*(t)  =  X(t)  -  A(t) 

2 
=  exp  (a   +  cut  +  a-t  )  -  exp(YQ  +  Yi"t) 


onto  the  event  stream  obtained  from  the  log-linear  intensity 
function. 

Once  it  is  determined  how  many  events  are  to  occur 
over  an  interval  with  intensity  function  A*(t),  i.e. 
N*  =  n',  it  is  necessary  tc  select  an  ordered  sample  of  n' 
variates  from  a  distribution  with  density  function 

X*(t) 


f*(t)  = 


[    u    X*(t)dt 


The  value  n'  will  be  small  compared  to  the  total  number  of 
events  in  the  original  non-homogeneous  Poisson  process. 
Therefore  the  need  for  realizing  efficiency  in  generating 
these  n1  ordered  variates  is  not  usually  as  crucial  to 
overall  algorithm  performance  as  is  the  need  for  efficiency 
in  producing  the  non-homogeneous  Poisson  process  from  the 
log-linear  intensity,  \ (t)  ,  in  Step  3.  However  if  the 
technique   for   obtaining   these   n*   ordered   variates    is 
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extremely  inefficient,  much  or  all  of  the  efficiency  gains 
realized  from  Step  3  above  could  be  lost  here.  (In  fact  an 
example  of  such  a  loss  of  previously  gained  efficiency  is 
documented  in  this  thesis;  see  Section  VI-B.) 

The  rejection  method  is  particulary  useful  for 
generating  random  variates  from  populations  with  continuous 
densities  that  are  bounded  and  which  are  concentrated  on  a 
finite  interval;  this  is  the  case  for  f*(t).  The  rejection 
method  and  an  algorithm  for  its  implementation  are  presented 
in  LEWIS  and  SHEDLER  [Ref.  11].  A  geometric  argument  will 
suffice  to  exhibit  the  principle  involved.  Consider  the 
density  function  f  (x)  ,  for  a  random  variable  X  on  the 
interval  (a,b)  in  Figure  13.  The  maximum  ordinate  of  this 
density  is  c.  If  another  function  g  (x)  =  c*  >  c  is 
constructed  then  the  density  function  f  (x)  is  enclosed 
within  the  rectangle  (a,0)  ,  (b,0) ,  (a  ,c*)  ,  (b,c*).  If  a 
point  within  this  rectangle  is  selected  at  random  it  will 
fall  either  under  the  density  function  or  above  it.  If  it 
is  under  the  density  curve  the  abscissa  of  that  point  is 
accepted  as  a  variate  from  the  population.  If  the  pcint 
lies  above  the  density  curve  that  point  is  rejected  and 
ancther  point  within  the  rectangle  is  selected  at  random. 
The  procedure  continues  until  n'  variates  have  been 
produced.  The  random  points  within  the  rectangle  are  easily 
produced  by  generating  two  independent  uniform  (0,1) 
variates  and  scaling  them   properly   resulting   in   a   point 

(x,y),  where  x  =  a  +  u-,(b  -  a)  and  y  =  U2#c*.  Then  if 
f(x)  >  y,  x  is  accepted  as  a  variate,  otherwise  it  is  not. 
After  n1  variates  have  been  obtained,  they  are  ordered  to 
yield  an  event  stream  for  a  process  with  intensity   function 

A*(t).  SDEEOOTINE  REJECT  employs  this  method  in  the 
program  for  Algorithm  B. 
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g(x)  =  c* 


*  X 


2  A-j  U 

The  density  f(x)   is  comD'etely  enclosed  by  the  rectangle 
(a,0),  (b,0),  (a,c*),  (b,c*). 

Procedure: 


1.  Generate  2  ^ndepende-r;  uniform  (0,1]  variates  u, 
and  lu. 

2.  Compute:  x  =  a  +  bu, ,  y  =  c*u?. 

3.  Plot  (x,y). 

4.  If  the  point  (x,y)  lies  under  the  density  f(x)  accept  x 
as  a  variate;  otherwise  go  to  Step  1. 

Example:  In  the  graphical  example  above  x,  would  be  accepted 

as  a  variate  whereas  x?  would  not  be  accepted. 

(Note:  Ideally,  c*  =  c  and  a  and  b  correspond 

to  the  lateral  limits  of  the  density  f(x).  This 

minimizes  ^ejection  reqion.) 

Figure    13    -      THE    REJECTION-ACCEPTANCE    METHOD    0?    VARIATE 
GENERATION    FROM    AN    ARBITRARY    DENSITY 
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Since  the  probability  that  a  point  on  the  abscissa 
will  not  te  accepted  as  a  variate  is  proportional  tc  the 
area  within  the  rectangle  that  is  not  under  the  density 
function,  it  is  advantageous  to  make  this  area  as  small  as 
possible.  Therefore  it  is  best  to  set  c*  =  c  and  set  the 
values  a  and  b  equal  to  the  end  points  of  the  interval  over 
which  the  density  function  occurs.  This  is  desirable  but 
not  necessary.  Figure  13  illustrates  the  more  general  case 
where  the  rectangle  is  larger  than  necessary.  One  may  be 
required  to  select  a  c*  >  c  if  c  is  difficult  to  determine. 
Seldom  would  a  and  b  not  coincide  with  the  lateral  limits  of 
the  density  however. 

It  is  obvious  from  Figure  13  and  the  preceding 
discussion  of  the  rejection  method  of  variate  generation 
that  it  is  the  "shape"  of  the  density  function  which  is 
critical  to  the  validity  of  the  method  and  not  the  fact  that 
it  integrates  to  one.  Therefore  any  scaled  density  would 
preserve  the  relative  shape  of  the  density  function.  As 
long  as  the  value  c*  is  at  least  as  great  as  the  maximum 
value  of  the  scaled  density,  the  rejection  method  will  yield 
valid  variates.  The  intensity  f unction  A  *  (t)  may  be  thought 
of  as  a  density  scaled  by  the  factor  y '  =  /  °A*(t)dt. 
Algorithm  B  uses  the  intensity  function  A*(t)  rather  than 
the  density  function  A*(t)/y'  in  the  rejection  routine. 
This  eliminates  a  division  step  for  every  function 
evaluation  required  by  the  method. 

5«   Merging  Event  Streams 

The  event  streams  from  the  intensity  functions  A  (t) 
and    A*  (t)   must   be  merged  to  produce  an  event  stream  from 
A  (t)  =   x(t)  +  x  *  (t)  .   In  the  present  case  there   are   two 
event   streams   (one   long   and  one  short)  which  are  already 
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ordered.  This  structure  allows  for  superposition  rather 
than  a  more  general  (and  more  time  consuming)  ordering 
procedure.  To  do  this  merging  job  SUBROUTINE  COLATE  is 
called. 


E.   SUMMARY 


This  section  has  presented  a  general  discussion  of  the 
more  significant  and  interesting  procedures  used  to 
implement  Algorithms  A  and  B  efficiently  in  FORTRAN.  The 
reader  is  referred  to  LEWIS  and  SCHEDLER  [Ref.  11]  for  a 
detailed  listing  of  steps  for  Algorithm  B.  Program  listings 
may  also  te  consulted.   These  may  be  found  after  Appendix  E. 
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V.   METHOD  OF  COMPARISON  OF  ALGORITHMS 


A.   GENEHAL 


Conclusions  drawn  from  the  results  of  comparing  the 
efficiencies  of  two  computer  programs  are  only  as  sound  as 
the  measures  of  effectiveness  upon  which  they  are  based. 
These  measures  of  effectiveness  are  always  somewhat 
arbitrary,  so  their  selection  should  be  based  upon 
compelling  reasoning. 

Even  after  reasonable  methods  of  effectiveness  have  been 
defined,  only  a  finite  number  of  trials  (usually  a  small 
finite  number)  for  a  given  set  of  circumstances  may  be  run. 
Therefore,  general  statements  such  as  "this  algorithm  is 
better  than  that  algorithm"  can  seldom  be  made  with  complete 
confidence.  Yet  if  the  incomplete  information  gathered 
reveals  certain  trends,  generalizations  hypothesized  from 
such  trends  may  warrant  a  high  degree  of  confidence. 

This  section  describes  how  the  two  algorithms  are 
compared  and  the  rationale  for  choosing  the  method  employed. 
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B.   MEASURES  OF  EFFECTIVENESS 


1 .   Computational  Speed 

Perhaps  the  most  popular  efficiency  measure  for 
computer  programs  is  their  speed  of  execution.  Speed  is  a 
natural  standard  as  long  as  computer  time  remains  a  costly, 
scarce  resource. 

Algorithms  A  and  B  were  compared  in  terms  cf  the 
mean  time  required  to  generate  event  streams  from  a  given 
set  of  intensity  functions.  It  should  be  noted  that  since 
the  number  of  events  in  each  event  stream  is  a  random 
variable,  the  time  required  to  generate  one  event  stream  is 
also  a  random  variable. 

Several  event  streams  with  different  intensity 
functions,  each  having  a  different  expected  number  of 
events,  are  produced.  Event  streams  for  each  of  these 
intensity  functions  are  replicated  many  times  by  each 
algorithm.  The  total  execution  time  required  for  all 
replications  is  the  "speed"  measure  of  effectiveness. 

The  number  of  replications  varies  with  each 
intensity  function.  For  example,  an  event  stream  with  an 
expected  number  of  events  less  than  200  was  replicated  100 
times  whereas  event  streams  with  over  1000  expected  events 
were  replicated  only  30  times.  In  both  cases  the  total 
number  cf  events  produced  was  of  the  same  order  of 
magnitude.  A  priori  it  seemed  reasonable  to  assume  that 
computer  time  required  to  execute  each  algorithm's  program 
would   be  rcughly  proportional  to  the  total  number  of  events 
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generated.  Therefore  in  designing  the  experiment  the 
product  of  the  expected  number  of  events  times  the  number  of 
replications  (E[Nt  ]  x  #  reps)  was  kept  roughly  constant  in 
order  to  keep  demands  on  computer  time  reasonably  under 
control.  (Note:  Program  execution  times  were  isolated  from 
compilation  and  linkage  times  for  the  purpose  of  the 
computational  speed  comparison.) 

2 .  Computer  Memory  Requirements 

A  seccnd  natural  efficiency  measure  for  computer 
programs  is  their  core  memory  requirement.  Both  programs 
were  written  with  the  goal  in  mind  of  conserving  as  much 
memory  as  possible  without  unduly  increasing  program 
complexity.  Core  requirements  reductions  could  most  likely 
be  accomplished  in  both  programs  through  more  sophisticated 
programming  methods.  Therefore  this  comparison  has  a  major 
weakness  as  a  measure  of  effectiveness.  Recognizing  this 
caveat,  memory  requirements  for  each  program  were   compared. 

3 .  Fidelity 

The  guestion  of  paramount  importance  is:  How  well 
does  the  program  simulate  the  true  process?  This  is  a 
difficult  question  to  answer  when  dealing  with  a  finite 
number  of  realizations  on  a  process  whose  theoretical 
properties  can  be  validated  only  after  an  infinite  number  of 
realizations. 

Several  methods  of  comparison  are  possible.  The 
simplest  is  that  cf  comparing  the  sample  means  and  sample 
variances  of  the  random  variable  whose  realizations  are  the 
number  of  events  generated  on  an  interval  by  a  given 
intensity  function.   This  random  variable  should  be   Pcisson 
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distributed  with  mean  equal  to  the  integral  of  the  intensity 
function  over  the  interval.  For  a  Poisson  random  variable, 
the  variance  equals  the  mean.  This  test,  although  useful 
for  determining  if  either  algorithm  varies  grossly  from 
expected  performance,  yields  no  information  about  hew  the 
algorithms  are  distributing  their  simulated  events  over  the 
interval.  It  is  necessary  to  look  at  discrete  segments  of 
the  interval  and  examine  the  distribution  of  the  number  of 
events  produced  in  each  segment  by  our  simulated 
ncn-homogenecus  Poisson  process. 

From  the  definition  of  a  non-homogeneous  Poisson 
process  (see  Section  II-B)  we  know  that  the  number  of  events 
observed  over  any  segment  of  the  interval  must  be  a  Poisson 
random  variable  with  parameter  equal  to  the  integral  of  the 
intensity  function  evaluated  over  the  interval.  3y  dividing 
into  several  segments  the  interval  over  which  the  Pcisson 
process  is  being  simulated,  the  number  of  events  in  each 
segment  can  be  observed.  Two  hypothesis  tests  may  then  be 
made  for  each  segment. 

The  first  test  is  the  dispersion  test  [Eef.  3, 
p.  158].  This  test  seeks  to  ascertain  if  the  number  of 
events  observed  over  each  segment  is  distributed  as  a 
Poisson   random   variable.    Let   n  ,   n  ,   ...f   n,   be   k 

-L      A  K 

observations  on  the  number  N   of  events  occurring  in  a  given 

P_ 

fixed   segment   p,   and   let  n  =  (n  +  ...+n, )/k.  If  N   is  a 

pi        k  p 

Poisson  distributed  random  variable,  then  the  statistic 

rk    ,      -  ,2 
)  •  ,  (n   .  -n 

d   =    1=1    P'1   P 

P         n 
P 

has  a  distribution  which  is  approximated  by  a  chi-squared 
distribution  with  k  -  1  degrees  of  freedom.  If  the  interval 
has  been  partitioned  into  m  segments,  then  by  simulating  the 
ncn-homogeneous  Poisson  process  k  times,  we   can   perform   ra 
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hypothesis   tests   at   a  selected  significance  level  a  .   The 

test  statistic  d    would   then   be   compared   to   the   (1-a) 

p 

percentile    of   the   chi-squared   distribution   with   k  -  1 

degrees  cf  freedom.   For  a  large  number  of  hypothesis  tests, 

say   1000,  we  would  expect  approximately  1000a  rejections  of 

the  null   hypothesis,   if   in   fact   the   N  's   are   Poisson 

P 
distributed. 


For  example  consider  the  non-homogeneous  Poisson 
process  over  the  interval  (0,100].  If  this  interval  were 
partitioned  into  100  unit  segments  the  number  of  events 
occurring  in  each  segment  should  be  Poisson  distributed.  If 
51  event  streams  are  simulated  over  the  interval,  a  value 
for  the  test  statistic  d  may  be  computed  for  each  segment, 
i.e. 

■51  ,      -  .  2 
d 


v51  , 

>j  ,  (n   . -n 

-     l      P'1   ? 


p         n 


Assuming   a  significance  level  a=  .05  and  comparing  each  d 

2  P 

to  the  critical  value  Xrn   =  67.5048,   we   would   count   the 

number   cf   test   statistics   such  that  d   >  67.5048.  If  the 

P 
number  of  events  in  the  segments  are  indeed   Poisson   random 

variables  we  would  expect,  on  the  average,  that  5  out  cf  the 

100  dp  values  would  exceed  the  critical  value. 

In  the  above  discussion  of  the  dispersion  test,  no 
mention  was  made  of  the  parameters  of  the  Poisson  random 
variables  (assuming  that  they  are  Poisson) .  This  is  because 
the  dispersion  test  for  Poisson  distributions  is  parameter 
independent.  Thus  it  is  necessary  to  perform  a  second 
hypothesis  test  to  determine  if  the  mean  number  of  events 
per  segment  is  equal  to  the  difference  of  the  integrated 
intensity  function  evaluated  at  the  right  and  left  end 
points  of  the  segment  respectively.  For  large  k  the  sample 
mean  n  is  normally  distributed.   Under  the  null   hypothesis, 
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a. 
i.€.   that  the  mean  of  ND  =/DX(t)dt  =  y^  (where  a.  and  a.  are 

F   aj_         P         i       3 
the  end  points  of  segment  p)  n  is  normally  distributed   with 

mean  yp  and  variance  Up/k-   The  test  statistic  is 

e  ■  |5P  "  V 
P     AT7K 

and  the  critical  value  is  the  (1-y  )  percentile  of  the 
normal  distribution,  za/2*  Again,  there  would  be  one 
hypothesis  test  for  each  segment  and  if  we  made  1000  such 
tests  we  would  expect  to  reject  the  cull  hypothesis 
approximately  1000a  times,  on  the  average. 

Continuing  the  example  above,  it  is  now  of  interest 
to  see  if  the  mean  number  of  events  on  each  segment  agrees 
with  the  theoretical  value  determined  by  the  integrated 
intensity  function.  This  time  assuming  a  significance  level 
of  a  =  .10,  the  critical  value  for  each  test  statistic  Cp 
would  be  z-ct/2  =  l  -645.  Each  c0  >  1.645  would  be  counted  and 
if  in  fact  tie  random  variables  have  the  hypothesized  means, 
we  would  expect  (on  the  average)  that  approximately  10  test 
statistics  cut  of  100  would  exceed  1.645. 

After  the  dispersion  test  and  hypothesis  test  for 
the  means  have  been  performed  on  event  streams  from  each 
algorithm,  the  test  results  may  be  compared  to  see  if  either 
algorithm  appears  to  have  null  hypothesis  rejection 
percentages  which  are  not  consistent  with  the  a  levels  of 
the  hypothesis  tests. 

(It  should  be  noted  that  during  the  development  of 
the  programs  histograms  of  the  event  streams  were  plotted  to 
determine  whether  or  not  they  "fit"  their  respective 
intensity  functions.  Both  programs  seemed  to  produce 
satisfactory  results  by  this  subjective  measure.) 
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C.   OTHEE  CONSIDERATIONS 


1 •   Intensity  Function  Category 

Algorithm  B  classifies  the  input  intensity  function 
into  one  of  six  categories  and  then  determines  what  action 
to  take  to  generate  the  appropriate  event  stream.  It  might 
be  expected  that  the  speed  of  Algorithm  B  will  be  affected 
not  only  by  the  total  number  of  events  to  be  generated  but 
also  by  the  form  of  the  input  intensity  function. 
Decomposition  of  the  intensity  function  into  four  sections 
will  prctably  require  more  time  than  if  it  must  be  divided 
into  just  two  sections. 

Algorithm  A  treats  all  intensity  functions  alike. 
Therefore,  it  is  important  to  compare  the  two  algorithms  for 
all  six  possible  categories  of  intensity  functions.  It  is 
possible  that  Algorithm  B  could  be  faster  than  Algorithm  A 
for  one  category  of  intensity  function  and  slower  than 
Algorithm  A  fcr  a  different  category  of  intensity  function. 
The  six  intensity  functions  selected  for  the  comparison 
represent  all  of  the  categories  defined  by  Algorithm  B,  and 
are,  in  fact,  those  intensity  functions  illustrated  in 
Figures  7  through  12  of  Section  III.  These  categories  are 
numbered  frcm  I  through  VI  and  correspond  directly  to  those 
listed  in  Bef.  11. 

2 •   Evaluation  of  c*  Value  in  Rejection  Routine 

The  importance  of  reducing  the  size  of  the  rectangle 
enclosing  the  intensity  function  A*  (t)  prior  to  beginning 
the  rejection  algorithm  was  mentioned  in  Section  IV-D.  If 
possible   c*  should  correspond  to  the  maximum  value  of  \*(t) 
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on  the  interval  (0,t  ].  For  Cases  I  and  II  (Figures  7  and  8) 
this   maximum   value  is  easily  determined  since  it  occurs  at 
one  of  the  end  points  of  the  interval.   For  Case  III  (Figure 
9)   the   maximum  values  of  \  *  (t)  and  \*{t)     both  occur  at  the 

J.  « 

point  b  fchere  the  interval  (0,tQ ]  has  been  partitioned  into 
two  subintervals  (0,b]  and  (bft  ].  For  Cases  IV,  V  and  VI, 
(Figures  10,  11  and  12)  the  maximum  value  of  A*(t)  is  not  so 
easily  determined.  LEWIS  and  SHEDLER  [Ref.  11,  p.  11]  state 
that  it  is  not  possible  to  find  this  maximum  explicitly. 
However  an  upper  bound  may  be  determined  from  the  values  k 
and  k  where  \(t)  <  k  and  \ (t)  >  k  on  the  interval  (0,tQ  ]  (or 
subintervals  (0,b]  and  (b,tQ  ]  for  Case  VI).  Then  by  setting 
c*  =  k  -  k  it  is  certain  that  X*(t)  <  c*  on  the  interval  or 
sutintervals  cf  interest.  Figure  14  which  illustrates  the 
rejection  routine  for  the  intensity  functions  for  Case  V  and 
VI  used  in  the  analysis,  reveals  that  the  c*  computed  in 
such  a  manner  may  not  be  very  efficient  as  an  upper  bound 
for  the  rejection  algorithm.  One  would  expect  some 
reduction  in  speed  efficiency  for  Algorithm  B  when  c* 
greatly  exceeds  the  maximum  value  of  X*(t) . 


3 •   Designated  Tolerance  Level 

Algorithm  A  requires  selection  of  a  tolerance  level 
which  determines  the  stopping  criterion  for  Newton-Raphson 
iterations.  The  speed  of  the  algorithm  could  be  expected  tc 
be  very  sensitive  to  this  tolerance  limit.  By  setting  it 
too  small  the  comparison  would  be  prejudiced  in  favor  of 
Algorithm  B.  By  not  setting  it  small  enough  the  fidelity  of 
the  resulting  event  stream  to  the  true  intensity  function 
would  be  impaired. 

It    was    determined   that   a   tolerance   level   of 
_7 
E[N.  ]  x  10     occasionally   demanded   that   the    computer 

discriminate   beyond   its   significant   digit   capability  in 
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single  precision.   This   resulted   in   the   failure   of   the 

Newton-Raphson   iterations   to   converge,  although  in  theory 

they  would   have   done   so   if   enough   precision   had   been 

— ft 
retained.    A   tolerance  level  of  E[N.  ]  x  10    seemed  to  be 

0 
unnecessarily  demanding   on   Algorithm   A   even   though   the 

computer   could  discriminate  at  this  level.   A  stopping  rule 

fixing   £  =  E[Nt  ]  x  10    was  the   compromise   that   insured 

reasonable   accuracy   of   the  simulated  event  stream  without 

unfairly  burdening  Algorithm  A. 
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(Acceptance  region  accounts  for  less  than 
4%  of  the  total  area  in  the  rectangle.) 
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Figure  14a  -  Graph  of  Rejection  Routine  for  the  Case  V 

Intensity  Function  Analyzed 
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Figure  14b    -    Graph  of  Rejection  Routines   for  the  Case  VI 

Intensity  Function  Analyzed 
Figure     14    -       ILLUSTRATION    OF    RS JECTION-ACCFPT ANCZ    REGIONS 
FOB    SAMPLE    CASS    V    AND    CASE    VI    INTENSITY    FUNCTIONS 
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VI.    RESULTS  AND  CONCLUSIONS 


A.    GENERAL 


Computer  programs  implementing  Algorithms  A  and  E  were 
compared  in  terms  of  efficiency  as  explained  in  Section  V. 
This  section  documents  the  results  of  the  comparison  and 
discusses  general  observations  of  the  author  concerning  the 
twc  competing  programs.  Also  recommendations  for  further 
study  concerning  Algorithm  B  are  offered. 


MEASURES  Of  EFFECTIVENESS  RESULTS 


1 .   S£eed 

Table  I  convincingly  illustrates  the  superiority  of 
Algorithm  B  over  Algorithm  A  when  computational  speed  is 
used  as  a  measure  of  effectiveness.  The  column  in  Table  I 
labeled  "Time  A/Time  B"  is  the  ratio  of  the  total  time 
required  ty  Algorithm  A  to  produce  a  specified  number  of 
replications  of  the  given  non-homogeneous  Poisson  process, 
to  the  time  required  by  Algorithm  3  to  produce  the  same 
number  of  replications  of  the  process.  For  all  six 
representative  intensity  functions  the  P cisscn-decoraposition 
and  gap  statistic  algorithm  is  at  least  twice  as  fast  as  the 
time-scale  transformation.  For  large  event  streams  the 
superiority  of  Algorithm  B  is  even  more  pronounced. 
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Table    I    -       COMPUTATION    TIMES    FOR    EVENT    STREAMS 
(IBM    System    360/67,     FORTRAN    IV    (G)     compiler) 
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The  most  obvious  disparity  revealed  in  Table  I  is 
the  difference  in  speed  efficiency  (of  Algorithm  E  with 
respect  to  Algorithm  A)  between  the  group  of  Cases  I,  II  and 
III;  and  the  group  of  Cases  IV,  V  and  VI.  In  the  former 
group  Algorithm  B  was  5  to  7.6  times  as  fast  as  Algorithm  A. 
For  the  latter  group  Algorithm  B  was  only  2.2  to  2.5  times 
as  fast  as  Algorithm  A. 

This  difference  between  the  two  groups  immediately 
made  suspect  the  value  c*  used  in  the  rejection  routine  of 
Algorithm  B.  For  the  first  group  c*  is  set  at  the  least 
upper  bound  for  \*  (t)  .  For  the  second  group  c*  is  not,  in 
general,  a  least  upper  bound  for  A*(t)  but  is  only  an  upper 
bound.  To  determine  whether  or  not  this  difference  in  the 
"quality"  of  the  c*  values  could  have  such  an  adverse  effect 
on  time  efficiency  for  Cases  IV,  V  and  VI,  these  cases  were 
simulated  a  second  time  with  new  c*  values  which  were 
empirically  evaluated  by  graphical  methods.  These  modified 
c*  values  were  set  close  to  the  least  upper  bounds  of  the 
respective  A* (t)  component  of  the  intensity  functions  for 
each  case.  Marked  time  efficiency  gains  were  observed  (see 
Discussicn  column,  Table  I)  indicating  that  significant 
efficiency  losses  can  be  sustained  due  to  the  method  of 
setting  c*  values  in  Cases  IV,  V  and  VI. 

Three  factors  appear  to  influence  the  degree  of 
speed  efficiency  gains  of  Algorithm  B  over  Algorithm  A. 
First  is  the  selection  of  the  c*  value  as  discussed  above. 
The  second  factor  is  the  percentage  of  total  events  which 
Algorithm  E  must  generate  from  the  rejection  routine.  For 
the  specific  intensity  functions  analyzed,  this  percentage 
ranged  from  U%  in  Case  V  to  2U%  in  Case  VI.  The  fewer 
events  required  from  the  rejection  algorithm  relative  to  the 
number  generated  from  the  gap  statistics  algorithm,  the 
greater  the  efficiency  of  Algorithm  3. 
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Finally  the  total  number  of  events  in  an  event 
stream  affects  relative  efficiencies  between  Algorithms  A 
and  B.  This  is  because  as  the  number  of  events  to  be  sorted 
by  Algorithm  A  grows,  the  less  efficient  its  sorting  routine 
becomes.  Sorting  requirements  for  Algorithm  B  are  greatly 
reduced  due  to  the  properties  of  the  gap  statistic  procedure 
used.  Therefore  as  total  events  increase,  so  does  the 
efficiency  advantage  of  Algorithm  B  over  Algorithm  A 
increase. 

2  .   Computer  M  emory.  Requirements 

Algorithm  A  requires  approximately  72,000  bytes  of 
core  storage.  Algorithm  B  demands  about  92,000  bytes.  Both 
programs  are  costly  in  terms  of  storage  due  to  the  number  of 
library  routines  used  and  the  need  to  store  up  to  5000  event 
times.  The  difference  of  20,000  bytes  would  not  be 
important  unless  computer  capacity  were  being  challenged. 
Memory  requirements  for  both  algorithms  could  be  reduced  by 
reducing  the  event  stream  capacity  of  the  programs.  A 
reduction  of  capacity  from  a  maximum  of  5000  events  to  a 
maximum  of  2000  events  in  each  program  would  result  in 
storage  requirements  of  54,000  bytes  and  68,000  bytes  for 
Algorithm  A  and  Algorithm  B  respectively. 

3 .   Fidelity 

Table  II  lists  the  results  of  six  series  of 
hypothesis  tests  of  the  types  discussed  in  Section  V-B. 
Three  series  of  dispersion  tests  and  three  series  of 
hypothesis  tests  for  the  mean  were  conducted.  Intensity 
functions  for  Cases  III,  IV  and  VI  were  selected  for  these 
tests. 
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Fcr  each  of  the  three  cases,  a  large  number  (500, 
900  or  1000)  of  each  type  of  hypothesis  test  was  performed 
for  a  fixed  level  of  significance  a.  The  number  of  times 
each  null  hypothesis  was  rejected  was  recorded  and  the 
statistic  a  =  {(number  of  re jections)/ (total  tests)}  was 
computed  for  the  series  of  dispersion  tests  and  the  series 
of  hypothesis  tests  for  the  mean.  If  the  null  hypothesis  is 
true,  then  a  is  an  estimate  of  the  significance  level  a. 

From  Table  II  it  appears  that  for  the  hypotheses 
tests  associated  with  Case  III  and  Case  IV,  a  gives  a  very 
good  estimate  for  a,  the  true  significance  level.  For  Case 
VI  some  disparities  are  observed.  For  Algorithm  B  the  a 
value  for  the  dispersion  test  (i.e.  hypothesis  test  that 
variates  are  Poisson  distributed)  is  0.039  whereas  the 
significance  level  is  a  =  0.01.  Algorithm  A  yields  a  much 
better  a  value  of  0.012.  For  the  hypothesis  tests  for  the 
mean,  both  Algorithm  A  and  3  yield  a  values  almost  twice 
that  of  a  • 

These  results,  though  interesting,  are  inconclusive. 
For  the  dispersion  test,  the  statistic 

lk    ,(n   --n  )2 

d  =  1=1  V'1     EL 

p        n 

P 

is  only  approximately  distributed  as   a   chi-sguared   random 

variable.    It   has   been   shown  by  FISHER  [Hef.  4]  that  for 

Poisson  random  variables  with  low   expectations,   hypothesis 

tests   based   on   the   chi-sguared   distribution  may  produce 

spurious   results.    Also,   we   would   expect   that   if   the 

distribution   of   d   is  only  approximated  by  the  chi-sguared 

p 
distribution,   that    the    relative    error    between    the 

approximate   and  true  distribution  is  greatest  in  the  tails, 

i.e.  at  high  percentile  values.   Case  VI  was  tested  at  an  a 

level   of   0.01.   It  is  also  the  intensity  function  with  the 
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fewest  expected  number  of  events  of  the  three  intensity 
functions  tested.  Partitioning  of  the  interval  into 
segments  produced  some  segments  which  had  a  low  expected 
number  of  events  (see  the  right  tail  of  the  density  function 
pictured  in  Figure  12).  The  results  for  Case  VI  rather  than 
indicating  any  serious  flaws  in  either  of  the  algorithms 
suggest  further  research  into  finding  better  ways  of 
hypothesis  testing  for  sparse  event  streams  and  highly 
discriminating  levels  of  significance. 

Eesults  of  the  hypotheses  tests  on  Cases  III  and  17 
indicate  that  both  algorithms  simulate  event  streams  that 
conform  to  the  desired  hon-homogeneous  Poisscn  processes. 
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Table    II    -      RESULTS    OF    HYPOTHESES    TESTS    FOR    FIDELITY 
(IBM    System    360/67,    FORTRAN    IV    (G)     compiler) 
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C.   GENERAL  OBSERVATIONS 


1 •    Programming  Ease 

The  programming  of  Algorithm  A  is  simpler  than  the 
programming  of  Algroithm  B.  Because  the  time-scale 
transformation  technique  handles  all  intensity  functions 
alike,  several  different  cases  need  not  be  identified. 
Algorithm  B  must  categorize  the  input  intensity  function 
into  one  of  six  categories  and  then  must  treat  each  category 
in  a  unique  way,  A  rough  indication  of  this  difference  in 
the  degree  of  complexity  is  the  number  of  instruction 
statements  required  by  each  program.  The  program 
implementing  Algorithm  A  required  142  instructions.  The 
program  fcr  Algorithm  B  required  364  instructions. 

2-    Exact  Method  Versus  Approximate  Method 

Algorithm  B  employs  an  exact  method  in  generating 
the  non-hcmogeneous  Poisson  event  stream.  It  is  exact  in 
the  sense  that  all  event  times  are  calculated  directly  and 
the  precision  of  those  calculations  are  limited  only  by  the 
physical  constraints  of  the  computer. 

Algorithm  A  employs  a  Newton-Raphson  iterative 
method  tc  approximate  event  times  on  the  interval.  The 
stepping  criterion  for  the  technique  is  limited  by  machine 
precision  considerations.  Also  the  stopping  rule  does  not 
give  a  firm  account  of  the  magnitude  of  error  that  can  be 
expected  in  the  event  times.  The  epsilcn  criterion  is 
applied  tc  the  value  of  the   integrated   intensity   function 
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and  not  to  the  abscissa,  (or  time  interval  axis) .  For  the 
integrated  intensity  function,  F(t),  it  is  conjectured  that 
given  a  function  value  F  (t  j_)  =  u-j_ ,  if  a  t  =  t»  can  be  found 
such  that  |F(t')  -  F(ti)|  <  e  that  t»  is  very  close  to  t  i- 
Although  this  is  a  reasonable  assumption  given  the 
well-behaved  nature  of  the  integrated  intensity  function, 
the  problem  cf  not  knowing  how  close  t*  is  to  t^  may  tend  to 
reduce  one's  confidence  in  the  algorithm. 

3 •   Initialization 

The  initialization  time  required  for  the  programs 
implementing  the  two  algorithms  has  not  yet  been  mentioned. 
Computation  time  comparisons  did  not  include  compilation  or 
linkage  times  required  for  the  two  algorithms.  These 
initialization  times  were  significantly  different,  being 
approximately  16  seconds  for  Algorithm  B  and  only  8  seconds 
for  Algorithm  A.  Should  simulation  of  only  one  or  two  event 
streams  be  required  it  is  likely  that  Algorithm  A  may 
actually  require  less  total  time  than  Algorithm  B.  The 
difference  between  the  initialization  times  could  probably 
be  reduced  if  the  programs  were  to  be  rewritten  in  Assembler 
language . 


D.   CONCLCSICN  AND  RECOMMENDATIONS 


1 .   Conclusion 

For  general  usage,  the  Poisson-decompositio n  and  gap 
statistic  method  of  LEWIS  and  SHSDLER  [Ref.  11]  is  clearly 
superior  to  the  time-scale  transformation  method  in 
generating   a   non-homogeneous   Poisson    process    with   a 
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degree-twc  exponential  polynomial  intensity  function. 

The  main  advantages  of  the  Poisson-decomposition  and 
gap  statistic  algorithm  are  its  speed  and  its  qualification 
as  an  exact  method  of  variate  generation. 

The  time-scale  transformation  method  enjoys  some 
advantage  in  core  storage  requirements  and  in  lower  program 
initialization  time.  These  advantages  are  not  sufficient  to 
overcome  the  relative  deficiencies  of  much  greater 
computation  time  and  uncertainty  about  the  accuracy  of  the 
approximating  mechanism  for  determining  event  times. 

2  •   Recommendations 

The  full  potential  of  the  Poisson-decomposition  and 
gap  statistic  algorithm  can  not  be  realized  until  it 
includes  a  better  method  for  selecting  an  upper  limit  value 
c*  for  Cases  IV,  V  and  VI  intensity  functions.  Algorithms 
for  choosing  a  c*  value  that  will  be  close  to  the  maximum 
value  for  the  function  \ *  (t)  should  be  developed  and 
incorporated  into  Algorithm  B.  A  single  variable  search 
technique  (such  as  a  Golden  Section  search  or  Bisection 
search)  for  finding  the  maximum  value  of  a  unimodal  function 
may  be  appropriate. 

The  computer   program   for   Algorithm  B   should   be 

rewritten   in   Assembler  language  in  order  to  reduce  program 

initialization  time.   The  program  could   then  be   developed 
into  a  library  routine  for  general  use. 

The  question  of  fidelity  of  the  simulated 
ncn-homogenecus  Poisson  process  to  the  true  theoretical 
process  should  be  investigated  further.  Of  special  inxerest 
is   the   question  of  how  well  the  algorithm  simulates  sparse 
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event  streams.  Methods  for  conducting  the  dispersion  test 
and  the  hypothesis  test  for  the  mean  at  very  high  levels  of 
significance  (i.e.  a  =  .01,  .005)  for  intervals  with  a  low 
mean  number  of  events  are  needed. 
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APPENDIX  A 

PROOF  OF  THE  VALIDITY  OF  SCALING  THE  INVERSE  DISTRIEOTION 

FUNCTION 


Proposition:      Let   T   be   a      continuous      random      variable      with 

distribution    function        F     (•)     ,    such   that 

T 


FT(t)     =    0 


m    A(t)     -    A(0) 
T  yQ 


0    <    t  <   fco 


FT(t 


t    >    t, 


Let  Y  =  U  F  (T)  such  that  T  €  (0,t  )  and   is  a   real   number. 
0  T  0 

Then   the   random  variable  Y  is  uniformly  distributed  en  the 

interval  (0, u  ) . 
0 

Proof:  F  (•)  maps  the  line  segment  [0,t  1  onto  the  interval 

[0r1].  If   F  (•)   is  strictly   increasing   on   (0,t0)   then 

F"1  (•)  exists   on   the   interval  (0r1)  and  maps  (0,1)  onto 

<0,t  ).  Now, 


Pv(y)  =  P[Y    y]  =  PU.F^T)    y] 


y0rT 


y 


=  P[F  (T)  <   -L-)] 

1      u0 


=  P[T  <  fI1  (-*-)] 


-  -T 


^0 


=  F  {F"1(-^-)> 

1  L       U0 


F  (y)  = 


_  Y 


0  <  y  < 


^0 
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Therefore, 


dFY(y) 

dy 


V*>  - 10 


0  <  y  <  u. 


and   I   is  uniformly  distributed  on  the  interval   (0,  u_) 
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APPENDIX  B 

ERROR  FUNCTION  AND  DAWSON* S  INTEGRAL  REPRESENTATION  CF  THE 
INTEGRATED  INTENSITY  FUNCTION   A  (t) 


The  standard  form  for  the  error  function  is: 

J      e     du   . 


/¥  0 


Dawson's  integral  is  usually  represented  as: 

-2   *    2 


t    /   e    du  • 
0 


Both  of  these  integrals  may  be  calculated  to  any  desired 
degree  cf  accuracy  by  using  the  exponential  series 
expansion, 


oo        n 
x  p       X 

6     =  Jn  nT 

n=0 


2  2 

u  ,         -u 

Thus,      e        and      e        may   be    written  as 


2  °°       ,    2,  n  °°        2n 

,U  p        (u    )  p       u_ 


*"      -      1 


n!  Ln    ni 

n=0  n=0 
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and 


2     "  ,   2,n     "  t    , ,n  -  2n. 

,~u   =   Y  (~u  )   _   V  (-D  (u   ) 

A    n!        n  n! 

n=0  n=0 


respectively.   Integrating  eu   yields 


t    2 


t   oo    2n 


/  eu  du  =  /   I     H 


du 


0   n=0 


t   2n 


I    /    a 

~  j.   n 


du 


n=0  0 


fc2n+l 

L  (2n+l)n! 
n=0 


-2 
Multiplyirg  by  t   results  in 


-2 


t   .,2 


0°    ^2n-l 


t  -  /   *u  *»  -  J0  ikmr 


and  the  series  representation  for  Dawson's  integral  is 
revealed.  This  series  will  converge  for  all  t,  although  the 
value  of  the  integral  increases  very  rapidly  as  t  increases. 
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Using  the  same  argument  it  is  easily  shown  that   the   series 

t   -u^ 
representation  for  fQ   e    du  is 

-  (.1>n  t2nH 

nio  <2n+1)n!      ' 


The  error  function  is  obtained  by  multiplying   the   integral 
by  the  constant  2//F   ,i.e., 

.        t  2  oo  n      2n+l 

2       r        -u       ,  2         v       (-D       t 

J       e  du   = 


r~     i  r-  n    (2n+l)n! 

/tt      0  /tt      n=0 


Although  the  series  expansion  may  be  used  to  calculate   both 

functions,    there    are   very   efficient  computer   library 

routines  available  for  computing  Dawson's  integral  and  the 
error  function. 

The  FORTRAN  supplied  procedures  ERF  and  DERF  [Ref.  15] 
evaluate  the  above  function  for  input  values  of  t. 

The  IMSL  (International  Mathematical  and  Statistical 
Libraries,  Inc.)  FUNCTION  MMDAtf  [Ref.  6]  evaluates  Dawson's 
integral  for  input  values  of  t. 


It    remains   to   be   shown  that   the   intensity   function 

2 

A  (t)  =  exp  (a,  +  a,t  +  at  }  may   be   integrated   over   the 

interval  (0,t  ]  using  either  ERF  (or  DERF)  or  HMDAW. 
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Consider, 

to        fco 

A(tn)-A(0)  =  /    X(t)dt  =  /    exp(a   +  a,  t  +  a  t  )dt 
u         0  0  z 


By   completing   the   sguare   of   the  exponent  the  expression 
becomes 


A(tQ)  -  /  °  expja0  -  jjAJ 


a. 


a2(t  +  2^}' 


(  A  (0)  =  0,   since   0   is  the  left  end  point  of  the  interval 
over  which  \  (t)  is  defined)  . 

For  a0>  0, 


a. 


A(tQ)  =  exp  jaQ 


4a. 


tn     (  a   )  2 

/  u  exp  {  •ST  (t  +  7^)\      dt 
0       (   Z 


2a. 


Letting, 


u  =  /a7  (t  + 


<*1. 


2  vv-    2a 


du  =  /a_   dt 


K '  =  exp  /  an  - 


a. 


0    4a. 


and   adjusting   the   limits   of   integration,  the  expression 
becomes: 


-£-   /   eu  du  = 


/o^T  a 


/   e    du  -  /   e    du 


/aT~   (  0 
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where, 


and 


a. 
a  =  (Zap 


a. 


b  -  </5I'     fco  +  2^r 


The  expression  above  can  be  rewritten  as  follows: 


i(t  >  - 


K 


a. 


f.  2.-2  rb   u2 

{b  b    J   e 

0 


2  "2  r 

du  -  a  a    J 

0 


2 

;U   du> 


Inside  the  brackets  are  two  Dawson  integral  expressions 
multiplied  by  constants.  Outside  the  brackets  is  just  one 
more  easily  determined  constant.  Using  MMDAW  twice  and 
multiplying  by  the  appropriate  constants  will  give  the 
desired  integral  value. 


For  a  <  0, 


a?  )   fc0 


A(tn)  =  exp  {    an  -  2-i  >  /    exp 


0    4a 


2  ;  0 


»2|(t  +  3^)2  i   dt 


a 


=  exp  <  a.  - 


0    4a 


Mr 

-  }    J         exp 


a-.   2 


t/y^T  (t  +  _i)]   dt 


Substituting  as  before  gives 


A(t, 
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/[  a 
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/I  a 


b   _  2 
/   e     du 
a 

/? 


-    b     2 
—  J   e    du 


/?    0 


/   e 


-u 


du 
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and   the   error   function   can   be   recognized   within    the 
brackets. 


Two   error   function   calculations,   a   subtraction   and   a 

multiplication   by   a   constant   are  sufficient  to  determine 

A(t  )  . 
o 
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APPENDIX  C 

COMPARISON  OF  GAP  STATISTICS  ALGORITHM  AND  TIME-SCALE 
TRANSFORMATION  ALGORITHM  FOR  SIMULATION  OF  NON-HOMOGENEOUS 
POISSON  PROCESSES  WITH  LOG-LINEAR  INTENSITY  FUNCTIONS 


LEWIS  and  SHEDLSR  [Ref.  13]  present  two  algorithms 
for  simulation  of  non-homogeneous  Poisson  processes  with 
intensity  function  X  (t)  =  exp  (g  +0  t)  .  The  first  algorithm, 
which  uses  a  time-scale  transformation  is  easy  to  employ 
since  the  inverse  of  the  integrated  intensity  function  is 
known  explicity.  The  second  algorithm  uses  gap  statistics 
to  generate  event  streams. 

Both  methods  are  exact.  That  isr  ether  than  the 
physical  precision  constraints  inherent  to  the  computer,  no 
approximations  are  necessary  in  generating  event  times. 
However  the  gap  statistics  algorithm  obviates  the  need  for 
ordering  an  array  of  variates.  Also,  the  gap  statistics 
algorithm  takes  advantage  of  a  fast  standard  exponential 
variate  generator,  (Random  Number  Generator  Package 
LLRANDOM,  Ref.  12),  thereby  avoiding  the  time  consuming 
computation  cf  taking  logarithms.  The  time-scale  algorithm 
must  perform  both  the  ordering  of  variates  and  the 
computation  cf  logarithms. 

Ecth  algorithms  were  programmed  and  tested  for 
computational  speed  efficiency.  The  gap  statistics 
algorithm  was  roughly  twice  as  fast  as  the  time-scale 
transformation  algorithm. 
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The  SUBROUTINE  NHPP2  uses  the  gap  statistics 
algorithm  within  the  Poisson-decomposition  and  gap  statistic 
algorithm  (Algorithm  B)  .  It  is  the  presence  of  SUBROUTINE 
NHPP2  which  markedly  reduces  the  ordering  requirements  of 
Algorithm  B  from  those  necessary  in  Algorithm  A. 
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APPENDIX  D 


POISSON  VARIATE  GENERATION 


A.   EACKGEOONC 

Acknowledgement:  The  content  of  this  appendix  is  a 
paraphrase  cf  portions  of  Chapter  11  of  an  unpublished  book 
by  J.    H.  AHEENS  and  U.  DIETER  [Ref.  1  ]. 


The    Pcisson   distribution   has   the   probability   mass 
function 


pi  = 


i! 


where  p.  is  the  probability  that  exactly  i  events  are 
observed  to  cccur  in  a  unit  time  interval,  given  that  events 
arrive  at  an  average  rate  \   . 


The  intervals  between  events  in  a  Poisson  process  of 
rate  \  are  independent  and  have  an  exponential  distritution 
with  mean  i/\,    i.e. 


f(x)  =  Xe 


-Xx 


The  protatility  integral  transform  method  can  be  used  to 
obtain  a  sample  of  exponential  variates  from  a  sample  of 
uniform  (0,1)  variates,  u  ,u  ,.  .  .  ,  since 


x,  =  -  y  ^n  ui  ' 


x. 


=  -  -  in   u2  , 


9a 


A  A  rate  event  stream  can  then   be   simulated   from   the 
exponential  variates  using  the  following  interpretation: 


1st  event  occurs  at  x 

1 


2nd  event  occurs  at  x   ♦  x 

1     2 

etc. 


This  property  yields  two  simple  methods  for  generating  a 
realization  en  the  Poisson  random  variable.  These  are  the 
additive  method  and  the  multiplicative  method. 


B.   THE  ADBI1IVE  METHOD  FOR  GENERATING  A  POISSON  7ARIATE 


The  number  k  of  units  arriving  in  a  unit  interval  in  a 
Poisson  process  with  rate  A  must  be  found.  This  will  be 
the  k  for  which 


and 


xn  +  x0  +•  •  •+  x,  <  1 
12        k  — 


Xl    +  X2  +'"+   Xk  +  Xk+1  >    1 


or  equivalectly,  for  which 


and 


-Jin    u,    -    Jin    u_    -•••-    Zn   u,     <_   A 


-In  u,    -   Jin   u2    -•••-    in   u,     -    in  u,     ,    >    X  (1) 


(Note:       Since    u.    <    1  Vi      ,    -In   u  .    >   0) 

Thus      by      using      uniform    variates,    computing   logarithms, 
adding    and    counting,    the   Poisson    variate   is      obtained.         The 
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problem  with  the  additive  method  is  that  if  \  is  very 
large,  say  2000,  it  requires  considerable  effort  to  generate 
a  single  Pcisson  variate.  Although  addition  is  extremely 
fast  on  a  computer,  the  computing  of  logarithms  is 
relatively  slow. 


C.   THE  MULTIPLICATIVE  METHOD  FOR  GENERATING  A  POI5SON 
VARIATE 


iultiplying  (1)  by  -1  gives: 


and 


£n(u,  *u2 u.)  _>  -A 


£n(ul*u2  uk'uk+l)  <  "X  ' 


or  equivalently 


and 


u ,  »u- u,  >  e 

12  k  — 


ii  •  n    •  •  •  •  •  n  •  ii     <  p 
ul  u2  k  uk+l   e 


So  by  generating  uniforms,  successive  multiplication  and 
counting,  the  Poisson  variate  k  is  produced.  The 
multiplicative  method  eliminates  the  need  for  computing 
logarithms.  However  trouble  is  encountered  for  large 
values  since  e  becomes  very  small  as  A  increases. 
Underflow  problems  occur  for  A  values  of  approximately  175 
and  larger  on  the  IBM  360/57  computer  system  using  single 
precision  word  lengths.  The  precision  problem  may  be 
overcome  by  partitioning  the  Poisson  random  variable  into 
the  sum  of  two  or  more  Poisson  random  variables.  However 
the  average  number  of  multiplications  and  uniform  variates 
required   is   still   proportional   to   A  .    More   efficient 
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methods  exist  such  as  the  Gamma  method. 


D.   THE  GAMMA  METHOD  FOR  GENERATING  A  POISSON  VARIATE 

(The  following  discussion  is  a  condensed  account  of  Ahrens 
and  Dieter's  presentation;  Ref.  1,  pp.  1 1-20,  21 ,  22) 

In   order  to   obtain   a   sample   k   from  the   Pcisson 
distribution   of   mean  y  ,   select   a   positive   integer   n 
(typically   n   is   a   little   smaller   than  u ) .  Then  take  a 
sample  x  from  the  Gamma  distribution  of  parameter  n. 

(1)  If  x  >  u  return  a 
sample  k  from  the  binomial 
distribution  with  parameters 
n-1  ,  u  /x 

(2)  If  x  <  y  take  a  sample 
j  from  the  Poisson 
distribution  of  mean  y-x  and 
return  k  •*-  n  +  j. 


The  sample  x  simulates  the  n-th  event  (arrival)  in  a 
Poisson  process  of  rate  1 .  If  x  >  u  then  there  are  n-1 
arrivals  in  the  interval  (0,x),  and  each  of  these  has  a 
probability  cf  y/x  of  being  below  u  [Case  (1) ].  If 
x  <  y  then  the  n  simulated  arrivals  are  all  before  y , 
and  the  sample  j  indicates  the  additional  events  between  x 
and   y  [Case  (2)  ]. 

(At  this  point  Ref.  1  includes  a   formal   proof   of   the 
procedure,  which  will  be  omitted.) 
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The  explicit  algorithm  relies  on  an  efficient  method  for 
sampling  frcm  the  Gamma  distribution.  (Such  a  method  has 
been  implemented  as  W.  R.  CHURCH  Computer  Center  library 
routine  GAMA  by  D.  W.  Robinson  and  P.  A.W.  Lewis  and  is 
documented  in  Reference  10.)  The  constant  n  is 
arbitrarily  chosen  as  n  =  [0.875  u  ]•  This  avoids  the  Case 
(1)  almost  completely  if  u  is  large  and  a  simple  method  for 
sampling  from  the  binomial  distribution  becomes 
sufficient:  the  Bernoulli  process  is  simulated  in  Steps 
8-10  below.  The  algorithm  for  sampling  from  the  Poisson 
distribution  with  mean  u  -  x  in  Case  (2)  is  the  simple 
multiplicative  method  explained  above  and  is  employed  in 
Steps  3-5.  However,  if  u  -  x  is  still  large  (larger  than 
the  "cut-cff  point"  c)  the  entire  procedure  is  iterated  with 
a  new  auxiliary  sample  x  from  a  Gamma  distribution  of  mean 
n»  =  [0. 875  (y  -  x)  ]  etc. 


Algorithm  -  The  Gamma  Methodx  Ahrens  and  Dieter 

1.  Initialize  k  ■*-  0 ,  w  «-  y. 

2.  If  w  >  c  go  to  6  (c  =  16  was  determined 
empirically  by  Ahrens  and  Dieter  to  be 
reasonable)  . 

3.  (Start   Case  (2))  .    Set  p  ■*-  1  and  calculate 

b  «-  e 


4.    Generate   u   and  set  p  *■  p«u.   If  p<b  deliver 
k. 
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5.  Increase  k  *-  k+1  and  go  to  4. 

6.  Set  n  «-  [dw]  where  d  =  7/8  (d  value  was  also 
selected  empirically) .  Take  a  sample  x  from  the 
Gamma  (n,1)  distribution.   If  x  >  w  go  to  8. 

7.  Set  k  =  k  +  n,  w  «-  w-x  and  go  to  2. 

8.  (Start  Case  (1))  .  Set  p  «-  w/x. 

9.  Generate  u.   If  u<  p  increase  k  «-  k  +  1. 

10.  Set  n  ♦■  n-1.   If  n  >  1  go  to  9. 

11.  Deliver  k. 

Ahrens  and  Dieter  claim  that  the  computation  time  for 
the  Gamma  method  grows  ultimately  like  a  +  B  In  y  • 
Computation  time  for  the  additive  and  multiplicative  methods 
grew  like  u  .  Therefore  for  the  large  u  values  typically 
demanded  in  the  simulation  of  non-homogeneous  Pcisson 
processes,  the  Gamma  method  is  signally  superior  in  terms  of 
computation  speed. 


99 


APPENDIX    E 


GRAPHICAL    PRESENTATION    OF    NEWTON-B APHSON    METHOD 


A.       GENERAL    PRINCIPLES 


(known) 


t.  (unknown)  t 

Figure  E-l 

Consider   the   function   F (t)   pictured  above  in  Figure  E-1. 

The  objective  is  to  find,   for   a   known   value   Uj_   on   the 

vertical   axis,   a   unique   corresponding  value  tj  such  that 

F(tj_)  =  Uj_ .   If  an  explicit  expression  for   the   inverse   of 

F(«)   exists,   the   calculation   may  be  made  directly,  i.e., 

t   =  F~^(u  ).   Otherwise,  numerical  methods  must  be  used  to 
i        i 
approximate  t  . 

i 
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G,(t)  =  F(t)   -  u. 


Figure  E-2 


Form   the    new    function      G . (t   )     =   F(t   )  -    u.       (see      Figure 

E-2).         In      effect,      the      horizontal   axis  has    been    displaced 

upward   a    distance    u..      Now    if   a    t*   can    be  found      such      that 

Gj  (t*)    =    0,      then    F(t*)    =    u.    and   t*   =    t.f  the   desired    value. 


Assume    that   an    initial    approximation    for    t*   can    be    made, 

say      tJ     .         If   G. (t*    )    and  G! (t!    )    =    g.  (t J    )    can   be    comDuted 
J        l  1  v   1    '  1111 

we    can    write 

G,(t!) 

»i(ti»  "  it;-  fen    ' 


where 


is  the  intercept  of  the  line  tangent  to  G.  (•)  at 


Gi  ^i  )  with  the  t-axis.   Solving  the  above   expression   for 

G  .  (t  • 


t'       gives: 
2 


t"    =    fc«    -      1  '  '1 
fc2        fcl     gTTE]T 


It    is    evident    from   the    graph    that    t'       is   closer    than    t' 
tc    the    unknown    t*. 
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G.(t)  =  F(t)  -  u. 
A 


Figure  E-3 


Repeating  the  procedure  (Figure  E-3)  ,  using  t£  as  a  new 
approximation  will  yield  t*  a  still  better  approximation  of 
t*.  In  general,  given  an  approximation  t£  has  been 
obtained,  a  closer  approximation  t£+,  can  be  calculated  by 
the  relationship 


fck+i  ■  fck- 


Gi'V 

g±(V 


It  can  be  shown  that  successive  approximations  will  converge 
to  the  value  t*  provided  G- (t)  satisfies  certain  criteria 
[Ref.  5,  p.  449,450], 

Since  G.  (t)  has  the  form  of  a  distribution  function  F (t) 
it  satisfies  all  the  conditions  for  convergence  except  for 
the  case  where  the  interval  [tl  , t . ]  contains  an  inflection 
point.  If  the  inflection  point(s)  can  be  identified,  this 
troublesome  situation  can  easily  be  avoided. by  handling  the 
function  G.  (•)  in  discrete  sections  over  the  discrete 
intervals,  (0,t»  ),  (t«  ,t»  ),  ...,  (y.ytJJ  )  where  the   t!j 
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(j  =  1,  .  .  .  ,  k)  are  values  such  that  Gj_(t".  )  are  all  the 
inflection  pcints  of  Gj_(«)  over  (0,tg].  (For  the  special 
case  under  consideration  at  most  one  inflection  point  will 
be  encountered.  It  is  the  maximum  or  minimum  of  the 
intensity  function.) 

Successive  approximations  are  computed  until 
IG^  (tj,  )  |  <  £  at  which  time  t»  is  assumed  to  be  close 
enough  to  t*. 


B.   INITIAL  AFPBOXIMATIONS 


The  problem  of  obtaining  a  good  initial  approximation 
for  each  t-  was  solved  in  the  following  manner  (see  Figure 
E-4  below)  : 


Figure  E-4 
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The   interval   (0 , to ]   was   divided   into   several   segments 
(0,x  h       (t  ,t    ]/   etc.   Then  the  values  F(i-:)  (j  =  1,2,...) 
were  computed.   Each  resulting  abscissa  and   ordinate   value 
was  stored  in  an  array,  i.e., 


t 

F(t) 

0 

0 

Tl 

Ht}) 

T2 

F(t2) 

• 

• 

fc0 

F(tQ) 

For  each  value  u   the  array  *as  searched  until  two   function 

values   were  located  such  that  F(x.)  <  u   <  F(x.,,).   Either 

t  .  or  t  .  .  -,  was  then  selected  as   the   initial   approximation 
J  3+1  rr 

for  the  value  t.  (i.e.  t-  , )  such  that  F(t.)  =  u-.  The 
function  G- (t  )  =  F(t  )  -  u-  was  then  formed  and  the 
Newton-Raphscn  iterations  performed  until  the  epsilon 
criterion  was  satisfied. 


104 


c  s 

c c 

c  9 

C              SUBROUTINE    MTDINV  C 

C  C 

C              PURPOSE  C 

C                          SIMULATES    A    NON-HOMOGENEOUS     POISSON    PROCESS    WITH  C 

C                          QUADRATIC    EXPONENTIAL    INTENSITY       FUNCTION    OVER    A  C 

C                          GIVEN    INTERVAL    USING    A     TIME-SCALE    TRANSFORMATION  C 

C                          OF     A    HOMOGENEOUS    POISSON    PROCESS.  C 

C  S 

C              USAGE  C 

C                         CALL    MTDINV    ( A, Al , A2, EL, ER, I S , 1 1 , M,  I ER )  C 

C  C 

C               DESCRIPTION    OF    PARAMETERS  C 

C                          A         -    CONSTANT    IN    INTENSITY     FUNCTION,  C 

C                          Al       -     1ST    DEGREE    COEFF     IN    INTENSITY    FUNCTION.  C 

C                         A2       -    2ND    DEGREE    COEFF     IN    INTENSITY    FUNCTION.  C 

C                          EL       -    LEFT    END    POINT    OF     INTERVAL.  C 

C                          ER       -    RIGHT    END    POINT    OF     INTERVAL.  C 

C                          IS       -    RANDOM    NUMBER    SEED.        ANY    INTEGER    WITH    MINE  C 

C                                         OR    LESS    DIGITS.  C 

C                          II       -    0    FOR    TIMES    OF    EVENTS.  C 

C                                           1    FOR    TIMES    BETWEEN    EVENTS.  C 

C                          M          -    RETURNEO    VALUE.    EQUALS    NUMBER    OF    EVENTS     IN  C 

C                                         NON-HOMOGENEOUS    POISSON    PROCESS.  C 

C                          IER    -    ERROR    COOE.     MAY    BE    WRITTEN    ON    DEMAND.        I ER  C 

C                                           EQUALS       THE    NUMBER       OF       TIMES       PROGRAM    WAS  C 

C                                           FORCED    TO    ABORT    SIMULATION    AND    START    AGAIN  C 

C                                          DUE    TO    EXCESSIVE    EVENTS    (MORE     THAN    5000).  C 

C  C 

C              COMMENTS  C 
C                          CALLING    PROGRAM    MUST    HAVE    A    COMMON    REGION,    BLOKltC 

C                          OF    DIMENSION     (5000).  C 

C  C 

C                            EXAMPLE:       DIMENSION    T(5000)  C 

C                                                        C0MM0N/BL0K1/T  C 

C  C 

C                          CALLING    PROGRAM    MUST    USE    THE    FOLLOWING    JCL    CARDS  C 

C  C 

C                                         //      EXEC       FORTCLG, IMSL=DP  C 

C                                         //FORT.SYSIN       DO       *  C 

c  c 

C                          TIMES    TO    EVENTS    OR    TIMES       BETWEEN    EVENTS    WILL    BE  C 

C                          STORED     IN    CELLS    T(l)     THROUGH    T(M).  C 

C  C 

c c 

c  c 
c 

SUBROUTINE  MTDINV  ( A, Al , A2 , EL , ER,  IS,  1 1 , M ,  I  E R) 

DIMENSION  TAU(5000),  BRKT(5000,2) 

COMMON  /BL0K1/  TAU/BLOK 2/ BRKT 

CALL  OVFLOW 

IER  =  -  1 
1  IER  =  IER+1 

BRKT( 1,1)  =0. 

BRKTC 1, 2)  =  0. 
C 

C      INTEGRATE  INTENSITY  FUNCTION 
C 

CALL  HELP  ( A, Al , A2 , EL , ER, P AR  ) 

PARI  =  PAR 
C 

C      GENERATE  POISSON  VARIATE  M 
C 

CALL  PVAR  ( IS, PARI, M) 

IF  (M.GT.5000)  GO  TO  1 
C 

C  SET    EPSILON    VALUE    FOR    NEWTON-RAPHSOM    STOPPING    RULE 

C 

EPS    =    FL0AT(M+1)*C. 00001 
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C  GENERATE    M    UNIFORM    VARIATES;    ORDER    AND    SCALE 

C 

CALL    RANDOM    <IS,TAU,M) 

CALL    PXSORT    (TAU,1,M) 
C 

DO    2    1=1, M 

TAU(I  )     =    TALK  I  )*PAR 

2  CONTINUE 
C 

C  DETERMINE     INTERVAL    PARTIONING    SCHEME 

C  FOR    SUBROUTINE    BNCHMK 

C 

12    =    M/4 

DIV  =  AMAXOCIO, 12 ) 

IDIV   =     INT(DIV+2) 

DELTA    =    (ER-ED/DIV 

ENDL    =    EL 
C 

C  CCMPUTE    ARRAY    OF    A8SC  ISSAS    AND    ORDINATES     FOR 

C  INTEGRATED    INTENSITY    FUNCTION 

C 

CALL    BNCHMK     < A, Al , A2 t ENDL, DELTA , IDIV ) 
C 

C  BEGIN    NEWTON-RAPHSON     ITERATIONS 

C 

K    =    0 

J    =    1 

3  K    =    K+l 

4  IF    (TAU(K).LT.BRKT( J  ,2)  )     GO    TO    5 
J    =    J  +  l 

GO    TO    4 

5  IF    ( ABS(BRKT( J,2)-TAU(K))  .LT,A6S(BRKT( J-L  ,2)-TAU( K)  )  ) 
*GC    TO    6 

T    =    BRKT1  J-1,1  ) 

ACD    =    ( BRKT(J-1,2 J-TAU(K) )/EVAL( A,A1,A2,T) 

GO    TO    7 

6  T    =    BRKT(  J,  1) 

ADD    =     < BRKT(J,2)-TAU(K) )/EVAL( A, A1,A2,T) 

7  IF    (ABS(AOD).LT.EPS)    GO    TO    8 
T    =    T-AOD 

CALL    HELP     ( A, Al, A2, EL,T,X) 

ADD    =    ( X-TAU(K) ) /EVAL( A, Al , A2,T) 

GO    TO    7 

8  TAU(K)     =    T 

IF    (K.EQ.M)     GD    TO    9 
GO    TO    3 

9  IF    (II.EQ.O)     RETURN 
C 

C  TIMES    BETWEEN    EVENTS    ARE    REQU  ESTED-CALCU  LATE    THEM 

C 

S    =    T  AU  (  1 ) 

TAU(l)     =    TAU(1)-EL 
C 

DO    10    1=2, M 

SI    =    TAU( I ) 

TAU(  I  )     =    TAU<  I)-S 

S    =    SI 
10    CONTINUE 
C 

RETURN 

END 
C 
C 

c 

C      SUBROUTINE  HELP  CALCULATES  THE  VALUE  OF  THE 

C      INTEGRATED  INTENSITY  FUNCTION  OVER  ANY  INTERVAL  (0,T) 

C      OVER  WHICH  THE  NHPP  OCCURS 


C 


SLBROUTINE    HELP     ( A, Al , A2, EL , ER, XX ) 
DOUBLE    PRECISION    MMDAWTBB,AA 
Z    =    SQRT(  ABSU2)  ) 
Y    =    (A1*Z) /<2.*A2) 
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AA    =    Z*EL+Y 

BB    =    Z*ER+Y 

CC    =    A-A1*A1/(4.*A2) 

CC    =    EXP(CC)/Z 

IF    (A2. LT.O.)    GO    TO    1 

Ql    =    DEXP( AA**2) *MMDAW( AA ) 

Q2    =    DEXP(BB**2 )*MMDAW(BB) 

XX   =    CC*(Q2-Q1) 

RETURN 
1    CC    =    CC*. 8862269 

XX  =  CC*(DERF(BB)-DERF< AA)  ) 

RETURN 

END 
C 
C 
C 

C      SUBROUTINE  PV<U  GENERATES  A  POISSON  (THETA) 
C      VARIATE,  M,  USING  THE  GAMMA  METHOD 
C 

SUBROUTINE  PVAR  (IS, THETA, M) 

DIMENSION  T(500O) 

COMMON    /BLOK1/    T 

K    =    0 

C    =    16. 

D    =    . 375 

1  IF     (THETA. LT.O     GO    TO    2 
GG    TO    5 

2  U    =    1. 

CTN    =    EXP(-THETA) 

MMAX    =    THETA+6.*SQRT(THETA) 

3  1=1 

CALL  SRAND  (IS,T, MMAX ) 
4-  U  =  U*TU) 

IF  (U.LT.CTN)  GO  TO  8 

I  =  1  +  1 

K  =  K+l 

IF  (  I.GT.MMAX)  GO  TO  3 

GO  TO  4 

5  NP  =  INT(D*THETA) 
AN  =  FLOAT(NP) 

CALL  GAMA  ( AN,IS,G,1) 

IF  (G.GT„THETA)  GO  TO  6 

K  =  K+NP 

THETA  =  THETA-G 

GO  TO  1 

6  U  =  THETA/G 
NP  =  NP-1 

CALL  SRAND  (IS,T,NP) 
C 

DO  7  1=1, NP 

IF  ITU  J.LT.U)  K  =  K+l 

7  CONTINUE 
C 

C      THE  VALUE  M  IS  ASSUMED  BY  THE  POISSON  (THETA)  VARIATE 
C 

8  M  =  K 
RETURN 
END 

C 
C 

c 

C  SUBROUTINE    BNCHMK    GENERATES    ARRAY    OF    ABSCISSAS 

C  AND    ORDINATES    OF    THE     INTEGRATED    RATE    FUNCTION. 

C  THIS    ARRAY     IS    USED    TO    DETERMINE     INITIAL     ESTIMATES 

C  FOR    NEWTON-RAPHSON    ITERATIONS 


C 


SUBROUTINE    BNCHMK    ( A, Al ,A2 , BINC , DELTA, ID  IV ) 

DOUBLE    PRECISION    AA,3B,MMDAW 
DIMENSION    BRKT(5000,2) 
COMMON    /BLOK2/    BRKT 
Z    =    SQRT(  /58SU2)  ) 
Y    =    ( A1*Z) /(2.*A2) 
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c 

c 
c 

AA    =    Y 

CC    =    A-A1*A1/(4.*A2) 

CC    =    EXP1CO/Z 

IF    (A2.LT.0.)    GO    TO    2 

Ql    =    OEXP( AA**2)*MM0AW( AA} 

DC    1     I=2,IDIV 
BINC    =    BINC+DELTA 
BB    =    Z*BINC+Y 

Q2    =    DEXP(BB**2)*MMDAW(BB) 
XX    =    CC*(Q2-Q1) 
BRKT( 1,1 )    =    BINC 
BRKTd  ,2)     =    XX 
1    CONTINUE 

GO    TO    3 
2    CC    =    CC*. 8862269 

DO    3    1  =  2,  IDIV 

BINC    =    BINC+DELTA 

BB    =    Z*BINC+Y 

XX    =    CC*(DERF<BB)-OERF< AA) ) 

BRKT( 1,1 )    =    BINC 

BRKT( I ,2)     =    XX 
3    CONTINUE 
C 

RETURN 

END 
C 
C 
C 

C      FUNCTION  EVAL  COMPUTES  THE  VALUE  OF  THE 
C      INTENSI TY  FUNCTION 


C 


FUNCTION  EVAL  (A,A1,A2,T) 

EVAL  =  A+A1*T+A2*T**2 

EVAL  =  EXP(  EVAL) 

RETURN 

END 
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c  c 

c c 

c  c 

C  SUBROUTINE    DEGTWO  C 

c  c 

C  PURPOSE  C 

C  SIMULATES    A    NCN-HOMOGEN EOUS     POISSON    PROCESS    WITH  C 

C  QUAORATIC    EXPONENTIAL    INTENSITY       FUNCTION    OVER    A  C 

C  GIVEN    INTERVAL      USING    THE       POISSON-0 ECOMPOSI TION  C 

C  AND    GAP    STATISTIC    ALGORITHM.  C 

C  C 

C  USAGE  C 

C  CALL    DEGTWO    (I S  ,  A  ,  Al , A2 , E L, ER , II , N,  I ER)  C 

C  C 

C  DESCRIPTION    OF    PARAMETERS  C 

C  IS       -    RANDOM    NUMBER    SEED.       ANY    INTEGER    WITH    NINE  C 

C  OR    LESS    DIGITS.  C 

C  A               CONSTANT    IN    INTENSITY     FUNCTION.  C 

C  Al       -    1ST    DEGREE    COEFF    IN    INTENSITY    FUNCTION.  C 

C  A2       -    2ND    DEGREE    COEFF    IN    INTENSITY    FUNCTION.  C 

C  EL       -    LEFT    END    POINT    OF     INTERVAL.  C 

C  ER       -    RIGHT    ENO    POINT    OF     INTERVAL.  C 

C  II       -    0    FOR    TIMES    OF    EVENTS.  C 

C  1    FD*    TIMES    BETWEEN    EVENTS.  C 

C  N          -    A    VECTOR      OF      LENGTH    5.    NU)     THROUGH      N(4)  C 

C  CONTAIN       NUMBERS       OF       EVENTS     FROM       VARIOUS  C 

C  COMPONENTS       OF       THE       DECOMPOSED       INTENSITY  C 

C  FUNCTION.       N<5)     CONTAINS       THE    TOTAL    NUMBER  C 

C  OF    EVENTS    IN    THE       NON-HOMOGENEOUS       POISSON  C 

C  PROCESS.  C 

c  c 

C  COMMENTS  C 

C  CALLING  PROGRAM  MUST  HAVE  A  COMMON  REGION,  HOLD,  C 

C  OF   DIMENSION  (5000),  AND  AN   INTEGER   ARRAY  OF  C 

C  DIMENSION  (5).  C 

c  c 

C  EXAMPLE:       DIMENSION    T(5000),    N(5)                                            C 

C  COMMON/HOLD/T                                                                       C 

C  C 

C  CALLING       PROGRAM          MUST       CONTAIN       THE       FOLLOWING    C 

C  ASSIGNMENT    STATEMENT:                                                                             C 

C  C 

C  M=N(5)                                                                                       C 

C  C 

C  CALLING    PROGRAM    MUST    USE    THE    FOLLOWING    JCL    CARDS    C 

C  C 

C  //      EXEC       FORTCLG,IMSL=DP                                                 C 

C  //FORT.SYSIN       DO       *                                                                    C 

C  C 

C  TIMES    TO    EVENTS    OR    TIMES       BETWEEN    EVENTS    WILL    BE    C 

C  STORED     IN    CELLS    T(l)     THROUGH    T(M).                                           C 

C  C 

c c 

c  c 

c 

SUBROUTINE    DEGTWO     ( I S , A , Al , A2 , E L , ER, II , N , I ER) 

DIMENSION    TIMES(5000),     T(5000),     N(5),    P<5) 

COMMON  /MIKE/  TIMES/HOLD/T 

CALL  OVFLOW 
C 

C      INITIALIZE  VARIABLES 
C 

P(l)  =  A 

P(2)  =  Al 

P(3 )  =  A2 

P(4)  =  0. 

P<5)  =  0. 
C 

DO  1  1=1,5 
1  N(  I)  =  0 
C 
C      IF  RATE  FUNCTION  IS  LESS  THAN  DEGREE  TWO, 
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C      USE  NHPP2  ROUTINE  ONLY 
C 

IF  (A2.EQ.O.)  GO  TO  2 

GG  TO  4 

2  GALL  NHPP2  (I  S, EL »ER, A, Al  t I  I »Nl  ,  IER) 
N(5)  =  Nl 

IF  (Nl.EQ.O)  RETURN 
C 

DO  3  I=1»N1 
TIMES* I)  =  T( I) 

3  CONTINUE 
RETURN 

C 

C      DETERMINE  COEFFICIENTS  FOR  MODIFIED 

C      DEGREE  ONE  RATE  FUNCTION 

C 

4  TEST   =    -A1/(2.*A2) 
TINT    =     ER-EL 

IF    (A1.GE.0.. AND. A2.GT.0.  )     GO    TO    5 
GO    TO    6 

5  B    =    A-A2*TINT**2 
Bi    =    A1+2.*A2*TINT 
GO    TO    10 

6  B    =    A 

IF    ((Al.LE.O. .AND. A2-LT.0.). OR. (Al.GT. 0-.AND.A2.LT. 0. 
*. AND. TEST. GE. TINT) )    GO    T3    7 
GO    TO    8 

7  Bl    =    A1+A2*TI,MT 
GO    TO    10 

8  IF    (A1.GT.0..AND. A2.LT.0. . AMD. TEST. LT . TINT)     GO    TO    9 
Bl    =    Al 

GO    TO    10 

9  Bl   »   Al/2. 
C 

C      MUST  THE  INTERVAL  BE  PARTITIONED? 
C 

10  IF  (Al*A2.LT.O.. AND. TEST. LT. TINT)  GO  TO  11 
ERNEW  =  ER 

GO  TO  12 

11  ERNEW    =    TEST+EL 
C 

C  GENERATE    DEGREE    ONE    NHPP    ON    INTERVAL 

C 

12  BB   =    B 

R  8  1     =    B  1 

CALL    NHPP2     <IS,EL,ERNEW,BB, BB1 i 0»N1, IER) 
N(  1)    =    Nl 

IF     (N(l J.EQ.O)     GO    TO    14 
C 

DO    13     1=1, Nl 
TIMES( I)     =    T( I) 

13  CONTINUE 
C 

C 

C  COMPUTE    LENGTH    3F     INTERVAL    AND    DETERMINE    VALUE 

C  OF    CSTAR    FOR   USE     IN    REJECTION    ROUTINE 


C 


14  Q    =    ERNEW-EL 
El    =    A 

E  2    -    A2*Q**2 

E3    =    Al *Q 

E4    =    A1**2/(4.*A2  ) 

E5    =    Al**2/< 2.*A2) 

IF    (A1.GE.0..AN0.A2.GT.0. )     GO    TO     15 

IF     (Al. LT. 0. . AND. A2.GT.0.. AND. TEST. GE. TINT)    GOTO     16 

IF    (Al.LT.O.. AND. A2.GT.0. .AND. TEST. LT. TINT)     GO    TO    17 

IF     (Al.LE.O.. AND. A2.LT.0. )     GO    TO     18 

IF     Ul.GT.O..  AND. A2.LT.0.. AND. TEST. GE. TINT)     GOTO    19 

CSTAR    =     EXP(E1-E4)-EXP(E1) 

GO    TO    20 

15  CSTAR    =    EXP(E1)-EXP(E1-E2) 
GG    TO    20 
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16  CSTAR  =  EXP(E1+E2+E3)-EXP(E1+E3) 
GO  TO  20 

17  CSTAR  =  EXP(El-E4-)-EXP(El-E5) 
GO  TO  20 

18  CSTAR  =  EXP( E1)-EXP(E1+E3+E2) 
GO  TO  2  0 

19  CSTAR  =  EXP(E1+E3+E2)-EXP(E1) 
C 

C      COMPUTE  INTEGRAL  OF  MODIFIED  DEGREE  TWO  RATE  FUNCTION 

C      OVER  INTERVAL 

C 

20  CALL    HELP     ( A. A  1 , A2 ,EL ,ERNE W ,PMTR ) 

PMTR    =    PMTR-(EXP ( B) *( EXP(  8 1*ERNEW )-EXP ( B 1*EL) ) ) /B 1 
C 

C  IDENTIFY    AS    FIRST    SUBINTERVAL 

C 

NOTE    =    1 
C 

C  GENERATE    REALIZATION    ON    POISSON     (PMTR)     VARIATE 

C 

21  CALL    PVAR     (IS, PMTR, M) 
IF    (NOTE.EQ.l )    GO    TO    22 
GO    TO    2  5 

C 

C      REJECTION  ROUTINE  USED  ON  FIRST  SUBINTERVAL 

C 

22  N(2)  =  M 
P<4)  =  B 
P(5)  =  81 

IF  (N(2)-EQ-0)  GO  TO  24 
CALL  REJECT  ( IS , EL ,C STAR, P , Q, N ( 2) ) 
C 

DO  23  1=1. M 

TIMEStN ( 1)+I)  =  T( I) 

23  CONTINUE 
C 

C      HAS  THE  INTERVAL  BEEN  PARTITIONED? 
C 

24  IF  (ERNEW.EQ.ER)  GO  TO  34- 
GO  TO  27 

C 

C      USE  REJECTION  ROUTINE  ON  SECOND  PART  OF  INTERVAL 

C 

25  N(4)  =  M 
P(4)  =  B 
P(5)  =  Bl 

C 

C  IF    NO    EVENTS    OCCURRED    BYPASS    REJECTION    ROUTINE 

C 

IF    (N(4).EQ.O)    GO    TO    35 

Q    =    ER-ELNEW 

CALL    REJECT    ( I S , ELNE W ,CST AR ,P , Q , N< 4) ) 
C 

C  COPY    TIMES     OF    EVENTS    INTO       'TIMES'     ARRAY 

C 

N4    =    N( 1)+N(2)+N( 3) 
C 

DO    26    1=1, M 

TIMES (N4+I)    =    T(  I  ) 

26  CONTINUE 
C 

C  GENERATION    OF    VARIATES    COMPLETE. 

C  GO    TO    ORDERING    ROUTINE 

C 

GO    TO    3  5 
C 

C  INTERVAL    PARTITION    WAS    REQUIRED.       MUST    NOW 

C  CONSIDER    SECOND    SUBINTERVAL 

C 

C  DETERMINE    COEFFICIENTS     FOR    MODIFIED 

C  DEGREE    ONE    RATE    FUNCTION 

C 
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27  IF     (A1.GT.0..AND. A2.LT.0. )     GO    TO    28 
B    =    A-A2*TINT**2 

Bl    =    A1+2.*A2*TINT 
GO    TO    2  9 

28  B    =    A+( A1/2.)*TINT 
Bl    =    Al/2.+A2*TINT 

29  ELNEW    =    ERNEW 
C 

C  GENERATE    DEGREE    ONE    NHPP    ON    INTERVAL 

C 

86    =    8 

B  B  1    =    B  1 

CALL  NHPP2  (IS, ELNEW, ER,BB,BB1,0,N3,IER) 

N(3)  =  N3 

IF  (N(3).EQ.O)  GO  TO  31 

N3    =    N< 1)+N(2) 
C 
C  TRANSFER    TIMES    BETWEEN    ARRAYS 


C 


DC    30    1=1, N3 
TIMES(N3+I)     =    T( I ) 
30    CONTINUE 


31  Q  =  TINT 
C 

C      DETERMINE  VALUE  OF  CSTAR  FOR  USE  IN 

C      THE  REJECTION  ROUTINE 

C 

E?  =  A?#Q**? 

E3  =  Af*Q 

IF  (Al.GT.O.. AND. A2.LT.0. )  GO  TO  32 

CSTAR  =  EXP(E1-E4)-EXP(E1-E5-E3-E2) 

GO  TO  3  3 

32  CSTAR  =  EXP(E1-E4)-EXP(E1+E3+E2 ) 
C 

C      CCMPUTE  INTEGRAL  OF  MODIFIED  DEGREE  TWO  RATE 

C      FUNCTION  OVER  SECOND  INTERVAAL 

C 

33  CALL  HELP  ( A, Al , A2 , ELNEW, ER ,PMTR ) 

PMTR  =  PMTR-(EXP ( 8 ) * ( EXP( B 1 *ER ) -EXP< Bl * ELNEW ) ) )/Bl 
C 

C      IDENTIFY  AS  SECONO  SUBINTERVAL 
C 

NOTE  =  2 

GO  TO  21 
C 

C      PARTITION  OF  INTERVAL  NOT  REQUIRED.   COMPUTE  TOTAL 
C      EVENTS  AND  SUPERPOSE  TWO  EVENT  STREAMS 
C 

34  N(5)  =  NQ)+N<2) 

IF  (N(Z).EQ.O)  GO  TO  38 

LBGN  =  N( 1)+1 

J8GN  =  1 

CALL  COLATE  ( LBGN  ,N( 5 )  ,  1 ) 

GO  TO  3  8 
C 

C      PARTITION  WAS  REQUIRED.  DETERMINE 
C      AMCUNT  OF  SORTING  NEEDED 
C 

35  N(5)  =  N( 1)+N(2)+N(3)+N(4) 

IF  (N(2 ).EQ.0.AND.N(4  ).EQ.O)  GO  TO  38 

IF  (N(4).EQ.O)  GO  TO  36 

IF  <N(2).EQ.O)  GO  TO  37 
C 

C      MLST  SUPERPOSE  FOUR  EVENT  STREAMS 
C 

LBGN  =  N(l )+l 

LFIN  =  N( 1)+N( 2) 

CALL    COLATE    ( LBGN , LFIN,  1 ) 

LBGN    =    LFIN+N(3)+1 

JBGN    =    LFIN+1 

CALL  COLATE  ( LBGN , N ( 5 ) , JBGN  ) 
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GO    TO    3  8 
C 

C  MIST    SUPERPOSE    FIRST    HALF    OF    ARRAY    ONLY 

C 

36  N2    =    N( 1)+N(2) 
LBGN   =    ISM  1)+1 

CALL    COLATE    (LBGN,N2T 1) 

GO    TO    3  8 
C 

C  MUST    SUPERPOSE    SECOND    HALF    OF    ARRAY    ONLY 

C 

37  KK    =    N(  1)+N(2)+1 

LBGN    =    N(l )+N(2)+N(3)+l 

LFIN    =    N(5) 

CALL  COLATE  (LBGN  ,N( 5  ) , KK ) 

GO  TO  3  8 
C 

C      ARE  TIMES  OF  EVENTS  OR  TIMES  BETWEEN  EVENTS  REQUESTED? 
C 

38  IF  (II.EQ.O)  RETURN 
C 

C  CALCULATE    TIMES    BETWEEN    EVENTS 

C 

S    =    TIMESd  ) 

TIMES(  1  )     =     TIMES(  l)-EL 

N5    =    N( 5) 
C 

DO    39    1=2, N5 

SI    =    TIMES(  I) 

TIMESd)     =    TIMES(  I  )-S 

S    =    SI 

39  CONTINUE 
RETURN 
END 

C 
C 

c 

C  SUBROUTINE    NHPP2     SIMULATES    A    NON-HOMOGENEOUS 

C  PCISSON    PROCESS    WITH    A    LOG-LINEAR    INTENSITY 

C  (RATE)     FUNCTION 

C 

SUBROUTINE    NHPP2     ( I S , EL , ER , A, Al ,  1 1  ,N , I ER  ) 

DIMENSION    T(500O) 

CCMMON    /HOLD/    T 


CALL    OVFLOW 
C 

C  INITIALIZE    VARIABLES 

C 

IER    =    0 

T  INT    =    ER-EL 

A    =    EXP(A+A1*EL) 
C 

C  IS    THE    POISSON    PROCESS    HOMOGENEOUS? 

C 

IF    (Al.EQ.O.)    GO    TO    3 

PAR    =     ( A*(EXP(TINT*A1  )-l. ) i/Al 

IF    (Al.GT.O.)     GO     TO     1 

IFLAG    =    3 

GO    TO    2 

1  A    =    A*EXP( TINT*A1 ) 
Al    =    -Al 

IFLAG    =    2 
C 

C  CCMPUTE    PARAMETERS    OF    BOTH    POISSON    RANDOM    VARIABLES 

C 

2  THETA    =    -A/Al 
GC    TO    W 

C 

C      COMPUTE  RATE  AMD  SCALING  PARAMETERS  FOR  HOMOGENEOUS 

C      PCISSON  PROCESS 

C 
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3  PAR    =    TINT*A 
I  FLAG    =    1 
AINVRS    =    l./A 

C 

C      COMPUTE  NUMBER  OF  EXPONENTIAL  VARIATES  REQUIRED 

C 

4  NMAX  =  PAR+6.*SQRT(PAR) 
C 

C      IS  THIS  A  HOMOGENEOUS  POISSON  PROCESS? 
C 

IF    (IFLAG.EQ.l)    GO    TO    17 
C 

C  GENERATE    REALIZATION    ON    POISSON     (THETA)     VARIATE 

C 

5  CONTINUE 

CALL    PVAR    { IS,THETA,M) 

IF    (M.EQ.O)    GO    TO    7 
C 

C  CALCULATE    TIMES    OF    EVENTS 

C 

CALL    SEXPON     (IS,T,NMAX) 

B    =    -Al 

V  =    0. 

JMAX    =    NMAX+1 
C 

DC    6    I  =  L,JMAX 
C 

C  HAVE    NUMBER    OF    EVENTS     EXCEEDED    THE    MAXIMUM    NUMBER 

C  THAT    THE    ARRAY    CAN    HOLD? 

C 

IF    (I.GT.NMAX)     GO    TO    8 

V  =  V+T(I  )  /UM-I  +  1  )*B) 
IF  (V.GT.TINT)  GO  TO  9 
T(I)    =    V 

IF    (I.EQ.M)    GO    TO    10 

6  CONTINUE 
C 

C  NO    EVENTS    OCCURRED 

C 

7  N    =    0 
RETURN 

C 

C      TOO  MANY  EVENTS  FOR  ARRAY.   INCREMENT  ERROR 

C      CODE  AND  TRY  AGAIN 

C 

8  IER  =  IER+1 
GO  TO  5 

C 

C  THE    NUM8ER    OF    EVENTS    OBSERVED    TO    OCCUR    IN    THIS 

C  NON-HOMOGENEOUS    POISSON    PROCESS     IS    'N1 

C 

9  N    =    1-1 

IF    (N.EQ.O)    RETURN 
GO    TO    11 

10  N    =    M 

11  CONTINUE 
C 

C      IS  THE  RATE  FUNCTION  INCREASING  OR  DECREASING? 
C 

IF  (IFLAG.EQ.3)  GO  TO  13 
C 

C  TIME    REVERSAL    TECHNIQUE     IS    NECESSARY 

C  DETERMINE    WHETHER    N    IS    EVEN    OR    ODD 

C 

ISIG   =    MOD(Nt2 ) 

NLOOP    =    N/2 

Nl    =    N+l 
C 

DO    12    1=1, NLOOP 

S    =   T(I ) 

T(I)     =    ER-T(N1-I ) 

T(N1-I)     =    ER-S 
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12  CONTINUE 
C 

IF    (ISIG.EQ.l)     T(NLOOP+l)=ER-T(NLOOP+l) 
C 

C  ARE    TIMES    OF    EVENTS    REQUESTED? 

C 

IF    (II.EQ.O)     RETURN 

GO    TO    15 

13  IF    ( II.NE.O)    GO    TO    15 
IF    (EL.EQ.O.)    RETURN 

C 

DO    14    1=1, N 
T(  I)    =    EL+T(I ) 

14  CONTINUE 
RETURN 

C 

C      CALCULATE  TIMES  BETWEEN  EVENTS 


C 


C 


c 
c 


c 


15  S  =  T(l  ) 


DO  16  1=2, N 
SI  =  T(  I) 
T(I)  =  T( I )-S 

S  =  SI 

16  CONTINUE 
RETURN 

C 

C      THE  POISSON  PROCESS  IS  HOMOGENEOUS 

C 

17  I  =  1 
U  =  0. 

CALL  SEXPON  (IS,T,NMAX) 

18  0=  U+T( I ) 

IF  (U.GT.PAR)  GO  TO  20 

I  =  1+1 

IF  (I.GT.NMAX)  GO  TO  19 

GO  TO  18 
C 

C      INCREMENT  ERROR  CODE 
C 

19  IER  =  IER+1 
C 

C  TRY    AGAIN    WITH    NEW    STRING    OF    VARIATES 


GO    TO    17 

20  N    =    1-1 

IF    (N.EQ.O)     RETURN 
IF    (II.EQ.l  )    GO    TO    22 

DO    21    1=1, N 

EL    =    EL+AINVRS*T( I) 

T(I )     =    EL 

21  CONTINUE 

RETURN 


22  DO    23    I =1 ,N 

T(  I)    =    T( I  )*AINVRS 

23  CONTINUE 
RETURN 
END 

C 
C 
C 
C 

c 

C      SLBROUTINE  PVAR  GENERATES  A  POISSON  (THETA) 
C      VARIATE,  M,  USING  THE  GAMMA  METHOD 


SUBROUTINE  PVAR  (IS, THETA, M) 
CIMENSICN  TC5000) 
CCVMON  /HOLD/  T 
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K    =    0 
C    =    16. 
D    =    .875 

1  IF    (THETA.LT.C)    GO    TO    2 
GO    TO    5 

2  U   *    1. 

CTN    =    EXP(-THETA) 

MMAX    =    THETA+6.*SQRT(THETA) 

3  1    =    1 

CALL  SRANO  (IS,T, MMAX) 

4  U  =  U*T(I) 

IF  (U.LT.CTN)  GO  TO  8 

I  =  1  +  1 

K  =  K+l 

IF  tl.GT.MMAX)  GO  TO  3 

GO  TO  4 

5  NP  =  INT(D*THETA) 
AN  =  FLOAT(NP) 

CALL  GAMA  < AN, IS,G,1 ) 

IF  (G.GT.THETA)  GO  TO  6 

K  =  K+NP 

THETA  =  THETA-G 

GO  TO  1 

6  U  =  THETA/G 
NP  =  NP-1 

CALL  SRAND  (IS,T,NP) 
C 

DO  7  1=1, NP 

IF  (T(  I  ).LT.U)  K  =  K+l 

7  CONTINUE 
C 

C 

C      THE  VALUE  M  IS  ASSUMED  BY  THE  POISSON  (THETA)  VARIATE 

C 

8  M    =    K 
RETURN 
END 

C 

C 

C 

C  SUBROUTINE    REJECT    GENERATES    AN    ORDERED    SAMPLE 

C  OP    GIVEN    SIZE    FROM    THE    SECOND    COMPONENT 

C  OF    THE    ORIGINAL    INTENSITY    FUNCTION 

C  USING    A    REJECTION-ACCEPTANCE    ALGORITHM 

C 

SUBROUTINE  REJECT  ( I S , EL ,CSTAR , P VEC ,Q , L ) 

DIMENSION  V(500),  PVEC(5) 

DIMENSION    T<5000) 

CCMMON    /HOLD/    T 

L20    =    L*10 

IF    (L20.GT.500)     L20=500 

LI    =    L+l 

K    =    1 

1  J    =    0 

CALL    SRAND     tIS,V,L20) 
C 

DO    2    1=1, L20 

J    =    J+l 

T(K)    =    Q*V<J)+EL 

J    =    J+l 

IF     (V( J) .LT.CALC( PVEC,T(K) )/CSTAR)    K=K+1 

IF    (K.EQ.L1)    GO    TO    3 

IF    U.GE.L20-1  )     GO    TO     1 

2  CONTINUE 
C 

IF     IK.LT.L)     GO    TO     1 

3  CALL    PXSORT    (T, 1 , L) 
RETURN 

ENO 
C 

c 
c 
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C  SUBROUTINE    HELP    EVALUATES    THE    INTEGRATED     INTENSITY 

C  FUNCTION    OVER    THE    INTERVAL     (EL,ER) 

C 

SUBROUTINE    HELP    ( A , Al , A2, EL »ER , XX > 

DCUBLE    PRECISION    MMDAW,BBtAA 

Z    =    SQRT( ABS(A2) ) 

Y    =    (AL*Z)/(2.*A2) 

AA    =    Z*EL+Y 

BB    =    Z*ER+Y 

CC    =    A-A1*A1/(4.*A2) 

CC    =    EXP(CC)/Z 

IF    (A2-LT.O.)    30    TO    1 

Ql    =    DEXP( AA**2 )*MMDAW( AA ) 

Q2    =    DEXP(BB**2)*MMDAW(BB) 

XX   =   CC*(Q2-Q1> 

RETURN 
1    CC    =   CC*. 8862269 

XX    =    CC*(DERF(8B)-DERF(AA) ) 

RETURN 

END 
C 
C 
C 

C      SUBROUTINE  COLATE  SUPERPOSES  TWO  ORDERED 
C      EVENT  STREAMS  OVER  THE  SAME  INTERVAL 
C 

SUBROUTINE  COLATE  ( LBGN , LFI NT JBGN) 

OIMENSION  TIMES(5000),  T(5000) 

CCMMON     /MIKE/    TIMES/HOLD/T 

I    =    JBGN 

J    =    I 

K   =    LBGN 

1  IF    (TIMES(  I  J.LT.TIMES(K)  )     GO    TO     2 
T(J)     =    TIMES(K) 

J    =    J  +  l 

K    =    K  +  l 

IF    (K.GT.LFIN)     GO    TO    3 

GO   TO    1 

2  TU)    =    TIMES(I) 
J    =    J+l 

I    =    1  +  1 

IF    ( I.EQ.LBGN)    GO    TO    5 

GO    TO    1 

3  11=    LBGN-1 
C 

DO    4   N=  1 1  1 1 
TU)    =    TIMES(N) 
J    =    J+l 

4  CONTINUE 
RETURN 

5  CONTINUE 

DC    6    N=KTLFIN 
T(N)    =    TIMES(N) 

6  CONTINUE 
C 

RETURN 

END 
C 
C 
C 

C  FUNCTION    CALC    EVALUATES    THE    SECOND   COMPONENT    OF 

C  THE    DECOMPOSED    INTENSITY    FUNCTION 

C  FOR    ANY     INPUT    VALUE. 


C 


FUNCTION    CALC    (P, ABSA) 

DIMENSION    P(5) 

X    =    P(l )+P(2)*ABS A+P(3 )*A3SA**2 

XX    =    P( 4)+P(5)*ABSA 

CALC    =     EXP( X)-EXP (XX) 

RETURN 

ENO 
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